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INTRODUCTION : 


The enclosed document is the final technical report on the 
work performed under NASA Cooperative agreement NCC2-158. 

This work covers the reports of three principal investigators 
who contributed to this work during the tenure of this project. 
This this report is a composite of three reports of: 

Dr. C.W. Bauschlicher , Dr. C.J. Nelin, and Dr. A. Komornicki 
The work where Dr. Komornicki was principal investigator was 
carried out by Dr. J. Almlof, of Oslo Norway. 



The following is the final report of Dr. C.W. Bauschlicher Jr. 
who was the principal investigator on this grant from September 
15, 1981 through June 30, 1982, at which point he resigned his 
position with the Institute. 

This report provided interim status on this grant, and was 
submitted to NASA technical personel in January 1983. 
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Introduction 


Recent research has focused on two main areas, 1) bonding in transi- 
tion metal oxides and 2) adsorption of CO and A1 and Ni. In both of 
these theoretical studies a major interest was to obtain a better under- 
standing of the nature of the bonding in transition metal containing systems. 
These calculations while useful themselves for the systems studied can 
also serve as a calibration of the methods employed when they are later 
applied to larger systems. These ab initio studies used self consistent 
field (SCF), multi-configuration self consistent field (MCSCF) and 
configuration interaction (Cl) methods in the treatment of the transition 
metal oxides and only the SCF method for the adsorption studies. In 
order to obtain reliable results, basis sets of reasonable quality were 
used. The contracted Gaussian basis sets used were of double zeta quality 
or better in the core and valence region or in the case of A1 in the 
valence region alone (the core having been contracted to single zeta 
using a segmented contraction). Further to be able to treat larger systems, 
for Ni the modified effective core potential (MEP) of Goddard (with 
the "4p" functions of Melius) was employed for the neighboring Ni atoms 
to the adsorption site and for the A1 neighbors used 3 electron ECP's. 

These studies have been carried out in collaboration with Paul S. Bagus 
(IBM, San Jose) and Charles W. Bauschlicher (Ames) . 

The transition metal oxide calculations focused on the two systems: 

CrO and MoO. These were part of a larger study (with Paul Bagus and Charles 
Bauschlicher) which also included NiO, PdO, CuO and AgO. Through this 
study both the effects of moving across the row (varying the number of 
d electrons) as well as moving down a column could be observed. In CuO 
there was no obvious direct involvement of the d electrons in the bond, 
in NiO there were excited states which had direct d involvement in the 
bond, whereas in CrO the lowest states formed Crd - Op bonds. These 
systems are primarily (Metal) + (0) with a charge of M).7 electrons. 

For Cr + the experimental d^ (ground state) - d^s^ separation is VL2, 000cm 
(SCF ^9,200cm ^) . Considering this large separation, it is surprising and 
still unresolved why many of the CrO states have as much d character as 


as they do. Often the states have d populations (near equilibrium) 

of approximately A. 3. For MoO the d populations for the high spin states 

5 4 1 + 

are approximately 4.6. (The experimental d -d s separation for Mo 
is ^12, 500cm ^ ; the SCF valve is ^18, 000cm \) At this point the changes 
between first and second row appear to be small, however the work on MoO 
is still quite incomplete as only high spin states have been calculated 
thus far. Unlike the dimers and Mo£ where Goodgame and Goddard found 
only an s-s bond for but a bond involving d for Mo 2 » CrO has clear 
d involvement in the bond. There are however larger s populations in 
many of the CrO states, especially the high spin states, than in MoO. 

One other striking feature of the CrO calculation is the radical bond 
shortening which occurs when configurations involving excitations from 

O 

02pn are included. At the MCSCF level the bond shortening was 0.1A 

O 

and at the Cl level it was M).2A. These configurations also contributed 
to bond shortening in NiO and CuO, though it was a much smaller effect 
in these systems. 

In the CO adsorption studies A1 and Ni were chosen as substrates so 
the effects of the valence sp electrons could be studied first alone 
(Al) and then in the presence of d electrons. A possible bonding mechanism 
involving mainly the valence electrons was first observed in a study of 
NiCO. It was seen that a state which had a valence configuration of 
4spO^ was very weakly bound whereas a state which had 4pfT^ was much 
more strongly bound. It was postulated that an increased p^ population 
increased bonding with C0(2tt ) whereas an increased spO population increased 
repulsion due to overlap with the CO 50. 

When clusters ranging in size from 1 metal atom to 9 metal atoms 
were calculated there was a changing balance between these two competing 
effects. With no neighbors (A1C0, NiCO) the state was either spO^ or 
pTT^ and unbound or bound (or barely unbound) accordingly. In the larger 
clusters the shape of the cluster tended to strongly affect the valence 
balance, i.e. lower layers reduced repulsion by "getting the o electrons 
out of the way", whereas neighbors on the same layer would compete for 
the electrons hence reducing the bonding. In all cases it appeared 



that the outer atoms donated charge to the central atom and bonding was 
dependent upon whether the o or n population increased. Although the 
role of the d electrons in the Ni clusters was never totally clear, when 
the clusters were treated as all 1 electron Ni atoms (the d’s now part 
of the MEP) there was no longer any CO backbonding, thus the d. electrons 
appear to certainly be involved in some subtle way. 

On the basis of the valenc population arguments it can be predicted 
that the order of binding with CO would be greatest for Na then decreasing 
for Mg and Al, i.e. the system with a smaller s population would be favored. 
Also transition metals would be favored as they have in general <1 valence 
sp electron. 

In the next two sections some detail on the oxide and adsorption 
calculation will be presented. 


Transition Metal Oxides: CrO and MoO 


The transition metal oxides are being studied because of their use 
as a basis for later treatments of chemisorption and oxidations of metal 
surfaces. Models for these processes will be substantially more involved 
than a "simple” diatomic. However the methods used must be calibrated 
for the diatomics and the features of the diatomic bonding must be under- 
stood before substantial surface related work can be undertaken. 

The contracted Gaussfan orbital basis set for Cr is that of Wachters 
including two function to represent the 4p orbital and also including 
Hay’s diffuse d function (14sllp6d/8s6p4d) . . The oxygen basis set is 
based on Dunning's double zeta contraction of Huzinaga's set, including 
also a diffuse p, a= 0.059 to represent the 0 character and ad a= 0.90 
(9s6pld/4s3pld) . . The Mo basis set is Huzinaga's set contracted using 
a generalized contraction (17sl3p9d/7s7p4d) .. 

For both Cr and Mo the ground states are ^S(d^s^) and the first excited 

states are ^D(d^s^) . For Cr + and Mo"*" the ground states are ^S(d^) and 

6 4 1 

the first excited states are D(d s ) . Both CrO and MoO and highly ionic 

molecules which are basically Metal"*" and 0 in character. For 0 in C 

1 4 2 3 

symmetry there are two possible states: o F (o hole) and on (n hole). 
When 0 is combined with either Cr + or Mo + , both ^’ 7 Z"*" and ^ ,7 II states 
can be formed. (The is the experimentally observed ground state). 

For both systems these four states are being calculated, with the 7 I + 
and 7 n states being treated at the SCF and Cl levels, and the and 
II states being treated at the MCSCF and Cl levels. The work on both 
systems is still ongoing. 

For the lower spin states two types of MCSCF calculations are considered. 
These are characterized by the configurations for which the orbitals are 
optimized. In the small calculation on the states 6 electrons are 

distributed among the 5 Cr + 3d orbitals and the 0 2pc orbital. In the 
large calculations there are 10 electrons and the 0 2pn orbitals are 
also included. The Cl performed was all singles and doubles from the 



six important configurations in the small MCSCF calculation. For both 
of the tt states broken symmetry calculations were performed, i.e. 

in C 2 V symmetry the two tt symmetries were not restricted to remain equi- 
valent. This was due to computation difficulties. However the separa- 
tion between the various states was much larger than the expected energy 
lowering due to symmetry breaking. The small MCSCF calculations for the 
H state had 6 electrons distributed among the 5 Cr3d orbitals and the 
0 2p7T x orbital, and the larger calculation had 8 electrons and the 2piry 
orbital is also included. 

For CrO at the MCSCF/SCF level the state ordering is ^11 (ground state), 
7 JI , ^£ + , 7 E + . At the Cl level the ordering of the upper states is the 
same and it is expected that the ^JI will remain the ground state (cal- 
culations are in progress). For the ^E + state r changes from 3.51a Q 

for the small MCSCF to 3.30a for the larger MCSCF and to 3.18a at 

o o 

the Cl level. The inclusion of excitations from the 02p7r orbital has 

a dramatic effect on the shortening of the bond length. At all levels 

of calculation near equilibrium the d population was ^.9. The net 

charge was ^0.8 and. the overlap was The two high spin states 

7 E + and 7 n both had considerably longer bond lengths ^r g = 3.5a o> which 

remained essentially unchanged from SCF to Cl levels. For both levels of 

theory the d population was approximately 4.2, and the charge was ^ 0.7 

at the SCF and ^4.3 and ^.6, respectively, for Cl. The (which is 

still uncompleted) has an r g 'V3.4a Q for the small MCSCF space, with a 

charge of M).6 and a d population of t 4.3. The tendency of these states 
4 

to have a large d character is surprising considering the reasonably 
5 4 1 

large d -d s separation, and is a question yet to be resolved. 



Adsorption of CO on A1 and Ni: Al x C0, Ni x CO 


To study the adsorption of CO on metals a cluster model approach was 
employed. With this method a system of atoms (cluster) is constructed 
which models many of the short range properties of the bulk system. 

While this technique gives reliable information on the equilibrium geometry 
of the adsorbate relative to the substrate, as well as the vibrational 
frequencies, information. on bond strengths is qualitative at best, however 
useful information on relative bond strengths can often be obtained. 

In these calculations the metal atoms were fixed at the bulk geometry 
and only the adsorbate-substrate distance was varied. In the present 
studies only one adsorbate molecule was present. Future studies hope 
to include several adsorbate atoms where the effects of lateral interactions 
between the adsorbates can be Investigated. These effects are important 
in understanding the differences between high and low adsorbate coverages. 

In these intial studies there were two main interests, 1) to understand 
the role of the s, p, and d electrons in the bonding and 2) to study 
the effect of varying cluster size and shape on the bonding (both in terms 
of nature and energy) . The cluster we designed to study the effects of 
having neighbors in the same layer only, having lower neighbors only, and 
then both effects competing simultaneously. In both systems only the 
(100) face and top sites are considered. 

Both A1 and Ni are face center cubic metals with lattice spacings 

0 O 

of 4.05A and 3.52A, respectively. For both systems 3 clusters were considered 
1) a single metal where there were no neighbor effects, 2) a five metal 
atom cluster which had one atom on the top layer and 4 lower neighbor 
atoms (1,4) [for Al a second 5 atom cluster was also considered which 

had 5 atoms in the top layer (5,0)], and 3) a nine atom metal cluster 

which had 5 atoms in. the top layer and 4 in the lower layer. For the Al 

clusters central atoms were all electron with a (4,3) CGTO basis and the 
neighbors used a 3 electron ECP with a (2,1) contracted basis. For the 
Ni cluster the central atom was all electron using Wachter's basis with 
p's and diffuse d (8s6p4d) . The neighbors used Goddard's one electron 


MEP with a (2,1) contracted basis. The C and 0 basis sets were (4,3) 

CGTO's for the A1 clusters and (4,2) CGTO's for the Ni clusters. 

3 1 

For the Ni clusters the results were NiCO A (with 4 spo ) unbound, 

Ni^CO (1,4) weakly bound and Ni^CO (5,4) unbound. In all clusters the CO 
it -backbonding (metal tt electrons being donated to the CO) was M).15. 
However only in Ni^CO was there substantial Ni 4pir character, and in keeping 
with the model of valence tt population aiding in the bonding, this is the 
only cluster which is bound. The cluster shape, 1 top 4 bottom, is 
consistent with the idea that the lower layer helps to reduce the sigma 
character by pulling it out of the way, whereas this cluster has no top 
neighbors to compete for the tt electrons and hence they can aid in bonding 
with the CO. 

For the A1 cluster only those clusters which had had a large 3piT 
character and a reasonably small a character were bound. These include 
states of A1C0, Al^CO (1,4) and Al^CO (5,4). For a given cluster differing 
electronic configurations could lead to the CO being either bound or un- 
bound depending upon the balance between the a and tt effects. 

Thus for both systems the role of the valence sp electrons appears 
to be understood, however although the d’s are certainly involved in the 
bonding of the CO to the Ni cluster their exact role is not yet clear. 

Finally for NiCO it has been observed that the inclusion of correla- 
tion has a considerable effect on the bonding, and it is hoped that these 
effects can be observed for larger clusters at a future date. 



The Chemisorption of Oxygen on Cu and N1 (100) 


The first phase of this work has involved the investigation of the 
diatomics NiO and CuO. The goal was to determine the extent of d involve- 
ment in these chemical bonds. This would determine the level of treat- 
ment required for the clusters. We found that there was little direct 
d involvement in the bonding, however their effect was to occupy space 
and indirectly influence the oxygen atom bonding by overlap effects. 

For Ni there was a little more direct d involvement than in Cu. In 
light. of the small d involvement, the d electrons were deleted from the 
Ni atom and small clusters of Ni atoms were used to model the surface. 

For the Cu atom an all-electron treatment was used as a further calibra- 
tion of the elimination of the electrons. A report of this work is 
summarized in the attached manuscript. 
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ABSTRACT: Model calculations indicate indicate two low-lying states for the interaction of 
O with Ni and Cu(100). The "oxide" state is «1 eV more stable and its R x for Ni(100) is 
consistent with the experimental value for Ni(100), indicating only one type of stable 
chemisorption bond. This agrees with XPS and EXAFS data for the observed LEED 
structures. Interpretations of EELS data in terms of a radical state for the p(2x2) and the 
oxide state for the c(2x2) structure must be seriously questioned. 

'Picscnt Addre>!>: N/-.SA Ames Research Center, Moffett Held, California 94035. 
**Mai!ing Address: 1 10. San Antonio Road, Suite 420. Mountain View, California 94043. 
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Upton and Goddard (UG), 1 using a 20 atom cluster to model the (100) surface of Ni, 
have shown that two low-lying electronic states exist for the Ni 2 gO cluster. One, termed the 
"oxide" state 2 has all oxygen orbitals doubly occupied either as lone pairs or as covalent 
Ni-0 bonds. The other, the "radical" state, has one singly occupied a orbital. On the basis 
of their results, UG and later Rahman, Black and Mills (RBM)^ accounted for the EELS 
spectra of p(2x2)0 and c(2x2)0 on Ni(100) in terms of the p(2x2) structures involving 
the "radical" state with the equilibrium distance, R , «0.9A and the c(2x2) structure 
involving the "oxide" state with ~0.26 A. In this communication we investigate 

theoretically the two states on model Ni(100) and Cu(100) surfaces. Our results and 
comparison to existing experimental data, with the exception of the EELS interpretation, lead 
us to conclude that there is only one stable chemisorption bond for O on Ni(100), the 
"oxide" state, with R a »0.89A. 

We model the (100) surface as a two layer, 5 atom cluster with 4 atoms in the top 

layer. The distance of the O atom above the 4-fold site is varied. For Ni we used the same 

modified effective potential (MEP) 4 as UG. Thus, only one Ni electron is explicitly included 

* 

in the calculation, the remaining 27, including the 3d shell, are present only in the one 
electron potential. The contracted Gaussian type orbitals (CGTO) representing the Ni 
contains 2s and lp functions. The p function is chosen to represent the Ni 4p orbital. The 
CGTO basis for O was (4s3p).^ It contains a diffuse p function and can represent both O 
and O'. For CujO. we used an all electron treatment. The O basis set is optimised for O' 
(4s3p) and the Cu is better than double ?eta in the 3d. 4s and 4p valence region. AJ1 basis 
sets are adequate to represent the electron di'tbbu‘ ; or.s rather well Symmetry and 
equivalence restricted self-consistent-ficld (SCF) methods were used to evaluate the 
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» 




wavefunctions. Correlation is included through the method of configuration interaction (Cl) 
and includes all single and double excitations from the SCF reference. 

The adequacy of the 5 metal atom clusters were tested by comparing to the N^qO 
results of UG. UG used only two s functions for Ni, so for comparison we deleted the 
Ni 4p functions. The "radical" state of Ni s O in C 4r symmetry is 4 B 2 ,...3aj4aJ le 4 2e 2 . The 
"oxide" state is 2 E,...3a^4a^le 4 2e 1 . The results of our SCF calculations at the equilibrium 
distance, R , of O above the Ni 4 square, the energy separation of the "oxide" and "radical" 
states, AE. and the calculated O vibrational frequency, a> e , are summarized in Table I. The 
results of UG, Walch and Goddard (WG)^ and the experimental EELS frequencies^ and 
LEED R^’s are also given. Agreement of our Ni^O calculations without the 4p functions 
with UG is very good. When, however, we put the 4p functions back into the calculations, 
R ± increases by 0. 1 1A (Table I) bringing it and the resultant <a e into agreement with the 
earlier calculation of WG on NijO. We conclude that the major difference between UG 
N^qO and WG NijO is the lack of the important 4p’s in UG and not a cluster size effect. 
This is consistent with our previous finding ** that the 4p’s contribute »50 percent of the 
cluster bond energy and our view that R x , u> t and the qualitative nature of the bonding all 
converge rapidly with cluster size. Note also that with the 4p's included, AE is only 0.1 eV. 

To test whether the two stales of oxygen are unique. to Ni, we invetigated Cu^O. Since 
an all electron calculation was performed, i! 3 lso serves as a check on the MEP 
approximation. We find the same two G states separated by a small AE and conclude that 
they are not unique to Ni and not an artifact of the MEP approximation. 

A Mulliken population analysis for our Ni s O (no 4p) cluster shows that the "radical' 
state has an oxygen charge of 0.67 and the "oxide" state has a charge of 1.32 in agreement 
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with UG. It is also well known that at the SCF level the computed electron affinity is 
-0.54 eV^, far from the experimental value of +1.462 eV, 1 * When correlation is included 
the electron affinity for O improves to +0.74 eV in our basis set. When we include 
correlation in our Ni s O calculation, we find that the energy of the "oxide" slate drops 
relative to the "radical" state by 0.79 eV to give a AE of 0.89 eV (Table I). This is 
anticipated because the large charge on the O atom in the oxide state should result in a larger 
correlation improvement. We also know that since we are computing only about half the 
experimental electron affinity of O at this level, further correlation improvements would 
lower the "oxide" state even further with respect to the "radical" state. 

One significant difference betweeen the Ni s O and Cu 5 0 oxide state results is the R x 
value, 0.74A for Ni and »0.98A (0.92A calculated SCF+ 0 . 06 A estimated for correlation) 
for Cu. For all-electron calculations on the diatomics NiO and CuO, we find 12 almost 

O 

identical bond-lengths with the computed value «0.1A longer than experiment (CuO). We 
also find that the MEP treatment for NiO gives a value ~0.2A shorter than the all-electron 
treatment. We therefore confidently expect that in the clusters the MEP calculated Ni-0 
bond length in the "oxide" state will be too short and the all-electron Cu-0 bond length too 
long. Averaging the two computed values (Ni and Cu; to give identical Ni-0 and Cu-O 

o o 

bond iengths. as is found for the diatomics, results in R^'s of 0.89A (Ni) and 0.82A (Cu). 

o 

The final conclusions from out calculations for NijO are, therefore, that R x is ~G.9uA lor 
the "oxide" state and it is >0.89 eV more stable than the radical state. 

UG suggested that the p(2x2)0 Ni(100) LEED structure represented formation of the 
"radical" O state on the surface because of agreement between their "radical" Rj and the 
LEFD determined value (O.hA) and because the calculated "radical" u' e was closer to the 
obsened FEES vibrational frequency 2 than that for the "oxide". Note that this choice is 
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made primarily because of the short R x calculated for the "oxide" by UG. We, however, 
find that the inclusion of the 4p functions, correlation, and an MEP bond length correction 
yield an "oxide" state with R x in agreement with experiment for the p(2x2)0 structure in 
addition to the ~1 eV greater stability of the "oxide" state. Our calculated’ « e values for 
the "oxide and "radical" states are practically identical when the 4p and Cl are included 
(Table I), though there is a large difference between the two without the 4p and Cl. 


UG attempted to simulate the higher oxygen covverage c(2x2)0 structure in their 
N^qO cluster by removing a charge from the cluster. This shortened the "oxide" state R ± 
even further (to 0.26A) and also gave an « e reduced relative to the 

"radical" state in agreement with the large reduction observed for the c(2x2)0 on Ni(100). 7 
UG therefore suggested that the higher coverage c(2x2)0 structure involved a drastic 
shortening of R ± and a switch to the "oxide" O state. RBM supported this opinion, since 
they wtre able to obtain a good fit to the experimental EELS vibrational frequencies using 
UG's u e values as input parameters for their lattice dynamical calculations. 


The interpretations of UG and RBM cast doubt on the LEED determined R x for the 
c(2x2)0 structured which is «0.9A, identical to the p(2x2)Od It is certainly true that the 
LEED data for c(2x2)0 and its fit to theory are poor and should be remeasured. There are 
however, several other experimental measurements which indicate that only a single type of 
O chemisorption state exists over the entire adsorption range up to oxide nuclealion (at 
~0.38 monolayers). 1 XPS shows only one 0( 1 s) B.F for the p(2x2)0 and c(2x2)0 
structures. Our orbital energies would imply a 1.6 eV B E difference between c>x;de and 
rauiCii> scales at '.hci; euuthbnum R . values, and still 1.2 eV ‘.or the oxide state at a R . as 
short nf 0 - V The ^ ork function increase ns a r :jnciion ccnerape chows only a smooth 
increase over the emirv chennsorpnon ranye, indicating the presence of only one type of 




O. 15 Finally, recent EXAFS data 1 *’ demonstrates both from adsorption edge characteristics 
and from EXAFS structure that the oxygen atoms in the p(2x2) and c(2x2) structures are 
identical in electronic character and R x . Their EXAFS derived value for R A 0.86 + 0. 07A is 
ia excellent agreement with our estimated value of 0.90A for the "oxide" state. 

We conclude that both the theory and experimental data for the Ni(100)/O system 
indicate that only one stable chemisorption state, the "oxide" state, exists and that R A for 
the state is «0.9A, as proposed in the original LEED data. We also predict that for 
Cu(100), oxygen will sit «0.83A above the surface. We feel that the Ni 2 oO + calculations 
of UG rather than simulating the c(2x2)0 structure is actually simulating oxide nucleation; 
hence the extremely short R A produced. The EELS data showing a large decrease in O 
vibrational frequency on going from p(2x2)0 to c(2x2)0 remains unexplained. Either there 
is some problem with the data, possibly involving some unrecognized oxide nucleation, or the 
claims^ that merely increasing O coverage without changing R A cannot, theoretically, result 
in the observed frequency decrease arc incorrect. We note that the experimental EELS 
changes on going from p(2x2)0 to c(2x2)0 are rather complex and continuous. 
Theoretically, we are directlydnvestigatir.g the effect of neighboring O’s on u t by placing 
several O’s on a large Ni cluster. Finally, with respect to the "radical state," we feel that it 
is a genuine low-lying state, but not the ground state at any coverage. It may have 
importance as an intermediate during reactions on Ni or Cu surfaces. 
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TABLE I 


Summary of Results 


Calculation/ 

Experiment 

Radical State 
R X (A) « e (cm -1 ) 

Oxide State 
R x (A) « c (cm _1 ) 

AE(eV) 

UG Ni 20 O a 

0.88 

371 

0.55 

218 

0.34 

RBM b 

• 

445 




UG Ni 20 O + a 

0.83 

355 

0.26 

266 

1.03 

REM b 


- 


295 


WG Ni 5 O c 



0.65 

290 


WG Ni 5 0 + c 

0.96 

300 




THIS WORK 

• 





Ni^O (no 4p) 

0.87 

374 

0.57 

251 

0.42 

with 4p 

0.98 

358 

0.68 

290 

0.10 

with 4p and Cl 

0.91 

317 

0.74 

320 

0.89 

Cu 5 0 SCF only 

!.37 d 

280 

0.92 

330 

0.21 

EXPERIMENT 






p(2x2)0 Ni(100) 



~0.90 e 

430 f 


c(2x2)0 Ni(100) 



"0.90 e 

3l0 f 



a Ref. I . 

b Ref. 3 , using the u; e in the line above as the single particle input parameter. 
c Ref. 6. The Ni^O calculation converged on the- "oxide'' state and we assume the Ni^O"*’ 
calculation converged to the "radical" state. 
c *Ave.r:;Ce of confige Titions used for the "radical" state 
e Rcf. 8, R determined from the theoretical fit to TV curves. 

*Rtf. 7. frequency of major peak, assumed to be surface-0 stretch. 
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Part I. Introduction 


For the period of July 1982 through January 1983, the research carried 
out under this grant has focused on the nature of metal adsorbate bonding. 
Two main areas were investigated in detail: 1) the adsorption of CO on 
Al, Na, Mg and Ni and 2) bonding in transition metal oxides, CuO, CrO and 
MoO. These theoretical studies have lead to a new understanding of CO- 
metal bonding as well as to a better understanding of the d involvement in 
the bonding in transition metal oxides. The research carried out, while 
useful in itself for the understanding gained, also serves as a calibration 
of the methods employed. This is necessary information if these methods 
are to be later applied to other systems. 

Part of this work has been accepted for publication and copies of 
the two papers are given in Appendices A and B. Other parts have been 
given as presentations at various meetings. This work was done in colla- 
boration with Paul S. Bagus (IBM Research Laboratory) and Charles W. 
Bauschlicher , Jr. (NASA Ames Research Center). 

This report is structured in three parts. Part II discusses the 
cluster model studies. Four areas are discussed in detail: hybridization 

and charging in cluster model studies, the adsorption of CO on metals 
(Al, Na, Mg), the adsorption of CO on Ni, and finally the use of Effective 
Core Potentials (ECP) in cluster model studies. In Part III the study 
of transition metal oxides, in particular CrO, MoO and CuO, is discussed. 



Part II. Cluster Model Studies of Metal-Metal and Metal-Adsorbate 


Interactions 


A) Introduction 

To study the nature of metal-metal and metal-adsorbate interactions 
a cluster model approach was employed. For this technique a bulk system 
is modelled by the use of a system of atoms (a cluster) which is construct- 
ed in such a way that it accurately predicts many of the local properties 
of the bulk system. While this technique gives reliable information about 
the equilibrium geometry of the adsorbate relative to the substrate, the 
vibrational frequencies, and the bonding mechanisms, detailed information 
on bond strengths is only qualitative in nature. Due to the finite size 
of the cluster, the metal atoms are held fixed at the bulk geometry and 
only the adsorbate-surface distance is allowed to vary. 

The primary emphasis of this study was to gain a better understanding 
of the mechanisms of CO-metal bonding. However we also analyzed the extent 
to which finite cluster effects (artifacts) can affect the calculation. 

Both the effects of changing the cluster size (number of atoms) and shape 
(inclusion of neighbors, next nearest neighbors, etc.) were investigated 
in detail. 

In part B of this section the nature of cluster effects as demonstrated 
through charging and hybridization will be discussed briefly. This work 
has been accepted by Chemical Physics Letters and the details are given 
in the attached paper (Appendix A). In parts C, D, and E the nature of 
CO-metal (A1 , Na, Mg), CO-Ni bonding, and the effects of ECP's will be 
discussed . 

B) Charging and Hybridization in the Finite Cluster Model 

For a cluster to accurately model the interaction of an adsorbate 
with a solid surface, the cluster wavef unction should have appropriate 
hybridization and polarization (i.e. close to that for the bulk material). 

To study the effect of cluster size on charging and hybridization 


clusters of Be, Mg and A1 were investigated. It had been thought previously 
that when all of the bonding site atoms have their nearest neighbors 
then r , w e and the nature of bond of an adsorbate to this cluster would 
be effectively converged. This rule of thumb appears to be true for the 
Mg and Be clusters, where for Be.^ and Mg^ the three types of atoms in 
the cluster, central, first layer edge and second layer are all essentially 
neutral. However for AI^q the central atom is negatively charged by 
0.63 electrons. When the size and shape of the A1 cluster was changed 
(Al^, Al^) there were further dramatic changes in the charging. It 
was also seen that the hybridization of the central atom is dependent on 
the cluster size and shape, eg. for Mg the hybridization as measured by 
the 3p population changes from 0.22 for Mg^ to 0.85 to Mg^. However in 
A1 there is similar central atom hybridization in Al^ , Al^ and Al^ 
cluster . 

Although it would be nice to have rules for the design of "good” 
clusters, it is now clear that the magnitude of finite cluster effects 
varies from material to material and even for clusters chosen to represent 
different faces of a crystal or different adsorption sites. Thus it appears 
necessary to investigate the question of appropriate cluster size separately 
for each chemisorption problem. More detail of this study is given in 
the attached paper (Appendix A) . 

C)The Bonding of CO to Metal Surfaces 

The commonly accepted idea of the bonding of CO in metal carbonyls 

is one of a donation and tt back-donation^ . The CO is oriented with C 

toward the metal atom and the CO 5o (C lone pair) orbital donates 

charge to and bonds with the metal. This a donation leads to a 

& 

back-donation into the unoccupied CO 2tt level which is also bonding. 

Angular dependent photoemission spectroscopy (PES) studies of CO on 

several transition metals (see Ref. 2 and references therein) have been 

interpreted to show that the molecular axis of the adsorbed CO is normal 

to the surface with C nearest the metal. Low energy electron diffraction 

(LEED) studies for C0/Ni(100) and C0/Cu(100) have further shown that CO 

3 

is in a head-on site directly on top of a surface atom . For CO adsorbed 



in this geometry the interpretation of the bonding, based on both 
4 2 

theoretical and experimental results , is the same as for the metal 
carbonyls, i.e. a donation and tt back-donation. In photoemission, ioni- 
zation of the CO 5 a derived orbital, for CO/Ni and on other surfaces, 

2 4-6 

shifts to a higher binding energy (BE) relative to the CO Itt level * ; 

these two levels of adsorbed CO have essentially the same BE. This has 

been taken as evidence that this CO a level is significantly, or even 
2 

primarily , involved in the bonding of CO to the surface. Theoretical 

studies have suggested that both the metal valence sp and the d electrons 

4 6 7 

are involved in the metal to CO bonding 5 * . 

In order to separate the effects due to the d and the sp electrons, 
we consider model studies of the interaction of CO with Na, Mg, and Al; 
they are all metals which do not have occupied d shells. In the next 
part the effects of d electrons will be discussed through the study of 
CO/Ni. Since we wish to apply our analysis to metals where the chemisorp- 
tion of CO has been studied rather extensively, we use the normal, head- 
on site geometry for CO as established for Ni, Cu, and other transition 
metals. The atoms treated differ primarily in the number of valence 
electrons but also in the symmetry character, 3s or 3p, of these electrons. 
The results are based on ab initio SCF cluster calculations for the 
interaction of CO with either one or four metal atoms; the clusters are 
denoted XCO and X, (1,3)C0 where X=Na, Mg, or Al. For each of the linear 
XCO clusters, we have considered two states: 1) a state with £ total 

symmetry where all the metal valence electrons have a symmetry, and 2) a 
state with II total symmetry where one of the metal valence electrons has 
tt symmetry. Thus we have one case, state (1), where there is no tt inter- 
action possible between the metal and CO and a second case, state (2), 
where this TT interaction is possible. By comparing the bonding in these 
two states we can more easily distinguish the nature of the interactions 
in a as opposed to tt symmetry . The X^(1,3)C0 clusters are chosen to 
simulate adsorption on the (1000) face of a hexagonal close packed (hep) 
structure for Na and Mg and the (111) face of face-centered cubic (fee) 

Al. The X^(l,3) clusters contain one first layer atom, the adsorption 
site, and its three nearest neighbors in the second layer of the crystal. 
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In these clusters there is a natural hybridization of the surface metal 

atom’s electrons into a and tt character relative to the CO axis; hence 

both types of interaction with the CO are possible. Polarization of the 

charge associated with the surface metal atom, either toward or away 

from the CO, is also possible and this may affect the nature of the 

interaction. For both the XCO and X.CO clusters the distance of CO 

4 . 

from the surface metal atom has been varied and the interaction has 

been studied as a function of this distance. 

Our analysis of the results of these cluster calculations leads 

to an interpretation of the bonding of CO which is quite different from 

the generally held view. We find that the interaction in the O space 

8 9 

is rather repulsive. A corresponding orbital analysis * shows that the 
orbitals of free CO are also, to a very large extent, occupied in the 
cluster, XCO or X^CO, wave funct ions . Thus there is very little, if any, 
a donation from CO to metal. We do, however, find that the metal charge 
in a symmetry hybridizes and polarizes away from CO so as to reduce the 
metal with CO 5a overlap and hence to reduce the repulsion. The interaction 

in the TT space is by contrast, quite attractive. The corresponding 

8 9 

orbital analysis * * shows that there is often a considerable transfer 

of charge from the metal TT electrons to CO and this is supported by the 

results of a more usual population analysis^. In the case of the X^(l,3) 

CO clusters, there is a competition for the bonding character of the tt 

electrons associated with the surface metal atom. These electrons (less 

than one) can bond this atom to the CO or they can bond it to the second 

layer metal atoms. In the later case the tt electrons contribute to the 

stability of the cluster (or the crystal surface which is being modelled) 

rather than to the stability of the adsorbed CO. The strength of the metal- 

CO bond can be related to the amount of tt charged transfered from metal 
* 

to CO 2tt . The net bonding of CO to the metal is the result of the sum 

of the a repulsion and the tt attraction. 

Thus we characterize the metal-CO bonding as a tt donation from metal 
* 

to CO 2tt coupled with a a polarization to reduce the repulsion in this 
space. Studies in progress indicate that this description of the bonding 
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applies for the interaction of CO with transition metals as well as for 

the simple metals discussed here. These studies will be reported else- 
11 12 

where * when they are completed. 

We have also examined BE shifts for the cluster levels derived from 

CO 5a, Itt, and 4a. We find that there is a large initial state shift 

of the 5a derived level to higher BE even for systems and metal-CO 

distances were the net interaction is quite repulsive. This shift is a 

consequence of the interaction and mixing of the metal and CO orbitals 

of a symmetry. However, it is not possible to relate the 5a derived BE 

shift to the strength of the chemical bond formed between the metal and 

2 

CO as has been attempted in the past . In particular we show that our 
results for A1C0 and Al^CO are consistent with the two peaked structure 
observed in the He I PES for low coverage of CO on Al(lll) by Schmeisser 
et al.^. Thus, even for the weakly bound adsorption of C0/A1 , M).21eV/ 

molecule, it may indeed be the case that CO adsorbs normal to the surface 
at a head-on site. 

D) Adsorption of CO on Ni(100) 

In the study of CO bonding to metals (Al, Na, Mg) we found that the 
dominant effect of the a electrons is repulsion whereas the bonding occurred 
in the tt space. We were then interested to see if this mechanism also 
applied to systems where d electrons were present. We hoped to find out 
what, if any, was the role of the d electrons in the bonding. Once again 
we are considering bonding to a head on site. The clusters studied were 
designed such that we could see the effects of having no substrate neighbors, 
NiCO, lower layer neighbors only, Ni 5 (l,4)CO, and finally both lower and 
top layer neighbors, Nig (5, 4) CO. It is expected that CO would bond more 
strongly to the Ni^ than to the Nig, cluster as in the Ni^ cluster the lower 
layer Ni atoms would tend to polarize the a electrons downward, hence 
reducing the repulsion and there would not be the competition for the 
tt electons between the CO and the neighboring surface atoms which would 
exist in the Nig cluster. There may also be different charging of the 
central atom (the adsorption) site which could effect the bonding. For 


the clusters the Ni atoms are fixed at bulk Ni lattice spacing, 3.52A, 

O 

the Ni-CO distance is 1.84A [the Ni-CO distance for Ni(CO) A ] and the C-0 
distance is 1.15A (the C-0 distance for free CO). In these clusters the 
adsorption site was treated as an all electron atom, and the neighboring 
atoms used 1 electron ECP ! s. (This approximation is discussed in more 
detail in Part E.) The basis sets used, contracted gaussian type orbitals 
(CGTO) , were [8,6,4] for the all electron Ni, [2,1] for the ECP Ni and 
[4,2] for the C and 0. 

NiCO 


Once again we were able to form two types of states, one with Ni 

9 13 9 13 

3d 4spa ( A) and the other with Ni 3d 4p7T ( $) . There is, however, a 

fundamental difference between these two states and those formed for the 


simple (non-transition) metal-CO (A1 , Na, Mg). In the J A state although there is 

* * 
no possible 4prr 2 tt bonding, it is now possible to have dTT 2tt bonding 

3 

occur. For the A state there is tt bonding, e.g. the CO tt population is 
0.11, however there is a comparable amount of repulsion. The net result 

is that the state is slightly unbound E = -0.09eV. On the other hand 

3 * * 

for the $ state, where both 4pTT 2 tt and dir 2tt bonding is possible, 

* 

we find C0(tt) = 0.68 and AE = leV. Thus while dTT ->• 2tt may contribute 

* 

to the bonding 4pir 2 tt is still the dominant bonding mechanism. 


Ni 5 (l,4)C0 

The Ni^ cluster is constructed with one Ni atom in the top layer (all 

electron) and four lower layer Ni atoms (ECP) . The state we consider has 
2 2 

a d Sx -y hole, however the different d holes in the a and 6 space are 

all quite close in energy. As expected this state is bound, AE = 0.50eV, 

3 

with C0(tt) =0.17. If we compare with NiCO A, which was unbound, we now 

have the opportunity to get bonding both from the dTT and the 4p7T (on NiCO 
3 

A only du could bond) . Furthermore comparing the populations of the Ni 

3 

in NiCO A with the central Ni in Ni^.CO, we see that the 4p7T population 
has changed from -1-0.01 (NiCO) to +0.31 (Ni^CO) , the pa population has 
remained similar +0.30 (NiCO) to +0.31 (Ni^CO) (however now this can be and is 
polarized out of the way) , and the s population has been considerably 
reduced from -0.03 (NiCO) to -0.56 (Ni^CO) . Thus the repulsion has been 
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considerably reduced while at the same time there are more tt electrons 
to become involved in the bonding. The net result being that the state 
is bound. Also, again, as in the simple metal-CO clusters, the valence 
sp electrons play a dominant role in the bonding and repulsion. Although 
the d electrons appear to be involved (this is best seen in NiCO) their 
role is by no means clear at the present time. 

Ni 9 CO 

In the Nig cluster four more ECP Ni atoms have been added as top 

layer neighbors, resulting in five atoms in the first layer and four in 

the second. As before the adsorption site is an all electron Ni atom. 

As predicted ^due to the competition for the tt electrons, they can either 

bond to the CO or to the neighboring Ni atoms. CO is unbound in this 

cluster, AE = -1.28eV. Although the CO tt population, 0.14, is approximately 

3 

the same as in NiCO A and Ni^CO, the Ni pTT population, -0.08, is more 
like that of NiCO than that of Ni^CO. 

As with the simple metal-CO we once again feel that the bonding of CO to Ni 
occurs in the tt space (predominately through the Ni valence 4TTp , however 

the role of the d’s is not yet established), and the a space dominately 
provides repulsion. 

E . ) The Use of ECP's in Cluster Calculations 

To treat larger clusters, especially those composed of heavy atoms, 
it becomes necessary to find a way to reduce the number of electrons which 
must be treated exactly. One technique rapidly gaining popularity is the 
use of Effective Core Potentials (ECP) . The philosophy behind this approx- 
imation is that in heavier atoms many of the electrons are inert with 
respect to the chemistry occurring, and as such they can be treated by a 
modelled potential which provides the appropriate interactions (mostly 
repulsive) with other electrons. 

As a part of the study on cluster effects we also examined the effects 
of the use of ECP’s both for the entire cluster and as neighboring atoms 
to the adsorption site. We looked at both A1 and Ni clusters . 

For the A1 cluster we focused on the A1^(1,3)C0 cluster, (111) face. 

The ECP used for this study was a 3 electron ECP in which the A1 3s and 3p 
electrons were left as valence electrons. The results of the SCF calculations 


are shown in Table I. Basically we find identical results between the all 
electron cluster and the cluster in which the adsorption site was all 
electron and the second layer neighbors were ECP’s. The small differences 
can possibly be explained as basis set effects (i.e. the use of the ECP 
introduces a slightly different basis set) . 

A second comparison was made for the Ni^(l,4)C0 cluster, (100) face. 

The ECP used in this case was a 1 electron modified ECP of Goddard. In 
this ECP only the 4s electron was left as a valence electron and furthermore 
the potential was modified to give the correct Ni atomic separations. 

For this study three different clusters were considered, one with 5 all 
electron Ni atoms, one with the adsorption site an all electron Ni atom 
and the second layer neighbors as ECP's, and finally one in which all 5 
Ni atoms, including the adsorption site, used ECP’s. The results are given 
in Table II. Once again the two clusters in which the adsorption site 
was an all electron Ni atom gave very similar results. However the all 
ECP cluster results were extremely different, and highly unphysical. 

[e.g. We have found that the bonding of CO to a metal or transition metal 
occurs in the tt space yet this cluster which is the most strongly 'bound 
of the three, D e = 1.33eV, has no tt bonding, C0(tt) = 0.03!] One problem 
with using a 1 electron ECP on the adsorption site is that the d electrons 
do not appear to be inert. It is probable that a cluster with a 10 elec- 
tron ECP (with the 3d and 4s electron as valence electrons) used for the 
adsorption site would lead to more reasonable results. 

From these two studies we feel that it is relatively safe to use ECP’s 
for neighboring substrate atoms, however great care must be taken if they 
are to be used for atoms directly bonded to the adsorbate. 



Part III. Study of Transition Metal Oxides 


A) Introduction 

The study of selected transition metal oxides was undertaken for two 
key reasons: 1) The results can be calibrated by existing experimental 

results, or in some cases they can provide a new interpretation of the 
experimental results; and 2) they can provide a foundation for future 
studies of oxygen adsorption and oxidation of metal surfaces. 

The transition metal oxides discussed in this section, CrO, MoO and 
CuO are part of a larger study (with C.W. Bauschlicher and P.S. Bagus) 
which also includes NiO, PdO, CuO and AgO. Through this study both the 
effect of moving across the row (varying the number of d electrons) as 
well as the effect of moving down the column (varying the spatial extent 
of the d electrons) could be analyzed. The results of this larger study 
are currently being organized and will be discussed elsewhere"^. This 
report will only discuss the specific results for CuO, CrO and MoO. (CuO 
will only be discussed briefly; a more detailed discussion of the results 
is given in the attached paper. Appendix B.) 

B) CuO - Predictions of a New Excited State 

Initially we choose to investigate CuO for two reasons: 1) because 

Cu is the simplest transition metal to treat computationally due to its 

15 2 2 + 

closed d shell, and 2) Huber and Herzberg had reported the X II— A E 
separation to be o/L6, 500cm \ which was considerably different from the 
E-IT separation, 350cm of the neighboring alkali oxide KO. As such 
CuO was viewed to be a prime candidate for d involvement in the bonding. 

We computed SCF and Cl wavefunctions for the ground and several low- 
lying states of CuO. One expects CuO to be a very ionic molecule with 
esentially Cu + and 0 character, and as such we classify the states accord 
ing to their free ion states. We find the ground state of CuO to be the 
X 2 n state with Cu + (3d^) , and we further find the first excited state 
2 E + [also Cu + (3d^)] to be located at ^8000cm a considerably different 
value than that given by Huber and Herzberg^. However this lower lying 




excited state has since been observed experimentally at 7865cm ^ ^ . 

We feel that the state reported in Huber and Herzberg was probably an excited 

E + state arising from Cu + in a 3d^4s* configuration. We predict that the 

+ 41 -1 

states formed from Cu (3d 4s ) will be 10,000 - 30,000cm above the ground 
state. 

2 2 + 

The computed separation of the X II and lowest E is now closer to the 

18 

separations obtained for the alkali oxides . The magnitude and direction 
of this separation can be explained by an analysis of the integrals between 
orbitals of the separated free ions. The details of this analysis and 
further details of the calculations are given in the attached paper 
(Appendix B) . 


C) CrO, MoO - A Clear Case of d Bonding 

In contrast to CuO where the lowest states had filled d shells, CrO 
and MoO can have the maximum number of open d shells. Once again the systems 
are expected to be highly ionic, essentially Metal + and 0 . From the 
separations given in Table III we would expect that the lowest lying states 
would have either Cr + or Mo + in a d~* configuration. However this simple 
model (looking at the ionic components) which worked well for CuO appears 
to break down when there are a large number of open d shells. There is however very 
little experimental information available for either CrO or MoO, hence 
most of the results obtained will be predictive in nature. 

20 

The contracted Gaussian orbital basis for Cr is that of Wachters 

including two functions to represent the 4p orbital and also including 
21 

Hay’s diffuse d function (14 , 11 ,6/ [8 ,6 ,4] ) . The oxygen basis set is 

22 

based on Dunning’s double zeta contraction of Huzinaga's set including 
also a diffuse p, a=0.059, to represent the 0~ character and a d, a=0.90, 
(9,6,1/[4,3,1]) . The Mo basis set is Huzinaga’s set contracted using 
a generalized contraction and also including a three term f function 
(17, 13, 9, 3/ [7, 7, 4,1]). 

14 

For 0 , in C symmetry there are two possible states: a tt (a hole) 


2 3 


aoy 


+ y 


and a tt (TThole) . When the 0 is combined with either Cr (3d ) or Mo (4d ) 


the resulting states are ^’^E + and 


Both of the quintet states have 
5vh 


the possibility of bond formation. For the "'E^ state the SCF configuration 
would be: 


(ID 


# 




o 




(da + 02po) 2 dTT 2 d6 2 02pTT 4 

and for the 5 JI state the SCF configuration would be 

(dn + 02pTT) 2 d7r02piT 2 da 1 d6 2 02pa 2 . 

For both CrO and MoO the high spin states, 2 II, are treated at the 

SCF level followed by a single reference singles and doubles CI(SDCI) . 

For the other states a MCSCF approach is used followed by a multi-reference 
SDCI. The MCSCF approach used is the Complete Active Space SCF (CASSCF) . 

In this approach the active electrons are allowed to form all possible 
spin and symmetry allowed distributions among the active orbitals. For 
both the and space two different active spaces were used. For the 

the active spaces in C^ v symmetry were: 

[3.1. 1.1] - 02pa, da, d5; di\; dTT ; d6 (6 electrons) 

[3. 2. 2.1] - 02pa, da, d6; 02p?r, dTT; 02p7T, dTT; dS (10 electrons) 
and for the II state the active spaces were: 

[2. 2. 1.1] - da, d5;02pTT, djr; d7T; d5 (6 electrons) 

[3. 2. 2.1] - 02pa, da, dS; 02p7T, d?r; 02p7T, dir; d6 (10 electrons) 
The resulting configurations from the small CAS calculations were bonding, 
anti-bonding and spin-recoupling. It turns out that the spin-recoupling 
terms are as, if not more, important than the anti-bonding terms. Thus 

for these systems a GVB wavefunction which neglects the recoupling terms 
would give very misleading results. For both the ^11 and the multi- 

reference SDCI is performed using the small CAS space as the reference. 

It should be pointed out that due to program limitations the ~*II 
states were done in broken symmetry (i. e. in C ^ symmetry the two tt symmetries 
were not restricted to remain equivalent) , and as such were not pure 
II states. However it is felt that the results from these broken symmetry 
states should well approximate those of the pure “*TI states. 

The quintet states, especially the for CrO, show very dramatic 
changes as the level of correlation increases (see Figure 1) . For this 

CrO state, on going from the small CAS to the SDCI , the equilibrium distance 

° -1 -1 
shortened by 0.32A, co changed from 420cm to 750cm , and D g increased 

by ^leV. One possible interpretation of these dramatic changes is that 
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for the small active space only the da-02pa bond was formed but then when 
the correlation was increased there was now also the possibility of some 
weak bonding in the tt space as well. This possibility is being investigated 
further. 

A summary of the CrO and MoO oxide results are given in Tables IV - 

VII and Figures 2-5. For CrO the state ordering at the CASSCF/SCF 

level is 5 * 7 II, 7 E + , ~*E + (Figure 2); for MoO the "*E + now drops below the 

7 n, 7 £ + to give the ordering ^n, ^E + , 7 II, 7 E + , (Figure 3). For both systems 

the Cl ordering is ^11, 7 II, 7 E + , (Figures 4 and 5). Thus we see for 

MoO the ^E + state was already reasonably well described by the small CAS 

calculations. For both CrO and MoO we find the ^II-^E + separation to be 

^6000cm ^ (see Tables IV and V), which is much the same as that for CuO, 

Or7500cm These separations have yet to be determined experimentally. 

For the CrO “*II state experimental information is available about r e 

and go and the Cl results agree extremely well (see Table V) . For both 
5 e + 5 

the E and II states of CrO and MoO the dominant bond formed is a d-02p 
bond (a or tt respectively) . 

A perplexing feature of these calculations is shown in Tables VI and 

VII. It had been predicted that the low-lying states of CrO and MoO should 

be formed from Cr + or Mo + in a d^ configuration, however many of the states 
4 1 

are more d s in character. A partial explanation for this effect can be 

5 4 1 

seen by examining the effect of a point charge on the d -d s separation 

of Cr + . If the energy of Cr + , in the field of a point charge at r = 3.252 a.u. 

4* 4 1 

is evaluated using the orbitals for a free Cr ion, we find that the d s 
2 2 -1 5 

state (dTT d6 ) drops 74cm below the d state. If the orbitals are allowed 
4 1 

to relax the d s state drops even lower. Certainly representing 0 as 

a point charge is an approximation and also the states are not pure ionic (e.g. 
4 1 

the d s type states of CrO and MoO have Q 'V/0.6) , however this effect 

4 1 

certainly is part of the explanation of the d s character. It is also 

5 4 1 

probable that in going from d to d s that the repulsion due to overlap 
of the Metal + orbitals with 0 orbitals is reduced. This possible effect 
is still under investigation. 

It should be stressed that the results presented here for CrO and MoO 
are preliminary in nature, and the computational data is still being 
analyzed and interpreted. 
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TABLE I 


Comparison between the population analyses obtained for an SCF 

calculation of an A1^(1,3)C0 cluster, (111) face, using four all electron 

A1 atoms and using one all electron A1 atom with three ECP A1 atoms (3 

electron). The calculations have C symmetry; the A1-C0 distance is 
o 3v 0 O 

1.98A, the A1 lattice spacing is 4.03A, and the C-0 distance is 1.15A. 

3 2 2 

The state compared is A 2 with valence occupation a e . Also included 
is a comparison of the results for an optimize A1-C0 distance. 




•> 


Top A1 Atom 
3s 

total 3p 3p 
r x r y 

3 >z 

Q 

Adsorbed CO 
a population 
tt population 

Q 


Properties for 
Optimize C0-A1 
(C0-A1) 


mm 


D 

e 


4 All Electron 


1 All Electron 


1.31 

1.40 

0.56 

-0.27 


1.28 

1.41 

0.51 

- 0.21 


-0.15 -0.15 

+0.40 +0.35 


-0.25 


- 0.20 


1.98A 


2.00A 


0.23eV 


0. 23eV 
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TABLE II 


# 

Comparison between the results for an SCF calculation of a Ni^(l,4)C0 
cluster, (100) face, using a) 5 all electron Ni atoms, b) 1 all electron 
Ni atom with 4 ECP Ni atoms (1 electron) , and c) 5 ECP Ni atoms (1 electron) . 
™ The Ni lattice spacing is 3.52A; the Ni-CO distance is held fixed at 1.84A, 

O 

and the C-0 distance is 1.15A. 



5 All Electron 

1 All Electron 4ECP 

5 ECP 

AE 

0.46eV 

0.54eV 

1.33eV 

co (a) 

-0.15 

-0.15 

-0.15 

CO (it) 

+0.19 

+0.16 

+0.03 


«> 
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TABLE 



Summary of calculations 
and Mo + ion. Separation, A, 

Cr 

for Cr atom 
is in cm ^ < 


E(3d 5 4s 1 ) 

E(3d 4 4s 2 ) 

A 

SCF 

-1043.2636 

-1043.2175 

10,120 

SDCI 

-1043.3221 

-1043.2866 

7,790 

EXP 3 


♦ 

8,090 




Mo 



E(4d 5 5s 1 ) 

E(4d 4 5s 2 ) 

A 

SCF 

-3975.4551 

-3975.3487 

23,350 

SDCI 

-3975.5268 

-3975.4414 

18,440 

EXP a 



11,830 


3 Ref. 18. 
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and Cr ion and for Mo atom 

nd the total energy is in Hartrees. 


E(3d 5 ) E(3d 4 4s 1 ) A 

-1043.0465 -1043.0046 9,200 

-1043.0950 -1043.0351 13,150 

12,280 

Mo^ 

E(4d 5 ) E(4d 4 5s 1 ) A 

-3975.2385 -3975.1557 18,170 

-3975.2948 -3975.2600 7,640 

12,800 


TABLE IV 


The CAS (MCSCF) and SDCI results for the lowest ^11 and states 
of CrO. The properties given are the equilibrium Cr-0 bond distance, r e> 
the vibrational frequency, U) e> and the dissociation energy D Also 
given are the energetic separation of the states, A and the correlation 
energy obtained for the SDCI wavef unction. These last two properties are 

o 

computed at R = 1.615A. 

^II A (cm -1 ) 

SCF r (A) 1.86 1.86 

e 

u (cm -1 ) 630 420 

e 

D (eV) 2.52 1.34 7210 

e v 


E (eV) 
corr 


1.615 

898.4 


See ref. 15. 
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TABLE V 


The CAS (MCSCF) and Cl results for the lowest "*TI and states 
of MoO. The properties given are the equilibrium Mo-0 bond distance, 
r^, the vibrational frequency, the dissociation energy, D^. Also 

given are the energetic separation of the states, A and the correlation 
energy abtained for the SDCI wavef unction . These last two properties 

o 

are computed at r = 1.77A. 


o 


SCF 

r e (A) 

to (cm ; 
e 

1.80 

560 

1.82 

730 



D (eV) 
e v 

2.13 

1.41 

6040 

SDCI 

r 

e 

1.71 

1.77 



03 

e 

1035 

830 



D 

e 

3.54 

2.63 



E (eV) 

corr 

7.29 

7.25 

6350 

£ 

Exp . 

r 

e 

— 

— 



0 ) 

e 

— 

— 



D 

5 .OeV 

— 



e 


a See ref. 15. 
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TABLE VI 


The resulting CAS and SDCI d population, the net charge on the Cr 


# 

atom and the net 
at r = 3.452 a.u, 

overlap population for 

5 n 5 i + 

the four 

7 n 

computed states 

V 

• 

CASSCF 

d 

4.26 

4.98 

4.25 

4.23 


Q 

0.63 

0.81 

0.64 

0.68 


overlap 

0.33 

0.30 

0.25 

— 

* 

SDCI 

d 

4.35 

4.89 

4.35 

4.32 


Q 

0.53 

0.78 

0.53 

0.58 


overlap 

0.46 

0.52 

0.32 

0.34 
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TABLE VII 


The resulting CAS and SDCI 4d population, the net charge on the Mo 
atom and the net overlap population for the four computed states of MoO. 




5 n 

5 Z + 

7 n 

V 

r (Mo-0) 

(a.u.) 

3.35 

3.35 

3.95 

3.95 

CASSCF 


d 

4.47 

4.89 

4.69 

4.80 


Q 

0.59 

0.82 

0.64 

0.74 


overlap 

0.42 

0.32 

0.24 

0.18 

SDCI 

d 

4.47 

4.90 

4.73 

4.84 


Q 

0.53 

0.73 

0.54 

0.64 


overlap 

0.54 

0.43 

0.30 

0.26 
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ABSTRACT: Cluster wavefunctions which have appropriate hybridization and polarization - 
lead to reasonable properties for the interaction of an adsorbate with a solid surface. 
However, for A1 clusters, we find that the atomic charge distribution is not uniform. 'The . ;. - 
finite cluster size leads to charges not representative for an extended system. This effect .... 
appears to be dependent on the particular material being studied; it does not occur in all 
cases. 
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The success of the cluster model to investigate some chemisorption problems relies on 
the local nature of the bonding of the adsorbate and metal surface atoms. Therefore, one 
must have a cluster large enough that all the atoms used to simulate an adsorption site on a 
metal surface have the same charge and hybridization and that these properties must be 
similar to that in the actual metal surface. We have investigated the question of the 
convergence of the cluster model for the chemisorption of H on a Be (1000) surface with 
respect to the cluster size. 1,2 We have considered clusters as large as Be^gf 14,8,14) (where 
the notation Be k (n,/,m) indicates n atoms in the first layer, l in the second and m in the 
third) and have concluded, that while the dissociation energy of the H from the surface, D e , 
converges very slowly with cluster size, properties such as the height above the surface, r e , 
the vibrational frequency normal to the surface, u g , and the nature of the bonding converge 
more rapidly. Our rule of thumb was that when all the atoms involved in the bonding have 
all their nearest neighbors, r e , « e and the nature of the bonding are effectively converged. 

We have shown that the number of neighbors is very important in determining the 
hybridization of the Be and Mg atoms,-* clearly a sufficient number of neighbors must be 
present so that the hybridization of the bonding atoms is similar to an atom in a metal 
surface. The cluster results support the view that the first nearest neighbors are needed for 
local properties to be converged. A simple example of this is H on Be (1000) at a three-fold 
site when there is no atom in the second layer below the incoming H atom; this is called the 
open site. Summarized in Table I are the r e , a> e and the atomic populations for the Be atoms 
directly involved in the H bonding at the open site for various size Be clusters. The results 
given are for minimal basis set SCF calculations. The smallest cluster, Be-j(3,0), includes 
only the directly bonding open site Be atoms. The largest cluster, Be^g, includes both first 
and second=nearest neighbors of these bonding Be atoms in the first three layers of the 
surface. The clusters and computational approach are described in detail elsewhere.*’ 2 For 



all but Be 2 2 and Be 36 , the bonding Be atoms are equivalent in the C 3v symmetry of the 
clusters. For Be 22 and Be 36 , the cluster has only C s symmetry and the average values of the 
populations for the nearly equivalent bonding atoms are given in Table I. The atomic 
valence electron population is decomposed into its 2s and 2p character. 

The variation of the computed « e is small among the various clusters. However, there 
are considerable changes in the r g . These changes can be related to the hybridization and, in 
particular, the decomposition of the p character into p z , normal to the surface plane, and p x 
or p y , in the plane. The directional character of the p orbital allows the formation of strong 
bonds and correspondingly short bond lengths. However, there is competition for use of the 
p z electrons to bond to the adsorbed H and to bond the first to the second layer of the 
cluster. 1 The H atom can also bond with the p x and p y Be electrons. However, this can 
only occur when the H atom approaches more closely to the surface. These two effects are 
the reason for the dramatic reduction in r e when a second plane is added to the cluster; note 
r g for Beg(6,0) and Be 4 (6,3). It is not surprising that one needs neighbors both in the 
surface plane and in the lower planes to induce the correct hybridization in the bonding 
atoms and to correctly represent the formation of metal-metal bonds. 

We have also investigated the chemisorption of oxygen on Al(lll) 4,5 For the bare 
Aljg cluster we found a charge of -0.29 electrons on the three central Al’s. 4 For this 
system, the oxygen-Al interaction is strong and O is negative because of a transfer of charge 
from the A1 atoms. Thus, the charge buildup at the central, bonding, atoms was assumed to 
only slightly effect the results; this is supported by the reasonably accurate results for r g . 
However, for weaker interactions, such charging of the bonding atoms might lead to serious 
errors. In this letter, we report on variations in charge distribution for several small clusters. 
Minimum basis sets were used throughout this work in large part because of the model nature 



of the calculations. The Be minimum basis is a 3 Gaussian type orbital (GTO) fit to Is, 2s 
and 2p Slater type orbitals (STO). This basis set is the same as that used previously and has 
been denoted Minimum basis lib. 1,2 Both the A1 and Mg basis sets are (633(s),42(p)) 
contractions of the Roos and Siegbahn (10s6p) primitive set 6 for A1 and their (10s4p) Mg 
basis with two p functions added to represent the 3p orbital.^ Symmetry and equivalence 
restricted SCF calculations were performed. The metal structure and metal-metal distances 
are taken from the bulk values.® The Molecule-Alchemy program system 9 was used 
throughout this work. 

Be and Mg form hep crystals and A1 has fee structure. For chemisorption at a directly 
overhead site on either the (111) surface of an fee crystal or the (1000) surface of an hep 
crystal, a 10 atom 2 layer cluster would appear to be appropriate. The central metal atom, 
the adsorption site, has its 6 nearest neighbors in the surface plane and its 3 nearest second 
layer neighbors. This Xjq( 7,3) cluster is shown in Figure 1. The total populations of the 
three distinct cluster atoms — central, first layer edge and second layer - are reported in 
Table II for Be 10 , Mg 10 and AIjq. The C 3v closed shell configurations used are: 

Be 10 : .... 7aj6e 4 la2 

Mg 10 : ....16a 2 19e 4 6a2 (1) 

Al 10 : .... 18a 2 20e 4 7a| 

For Be 10 and Mg 10 , the charges on the three types of atoms are essentially the same and 
essentially neutral. This is not true for Al 10 where the central atom is negatively charged by 
0.63 electrons. The gross atomic populations may be misleading because of the arbitrary 
assignment of overlap populations to atoms. However, the different behavior for A1 from Mg 
and Be stroTfgly suggests that the central A1 is indeed negatively charged. In order to further 
test the population analysis results, we examined the Is orbital energies; these are also shown 


in Table II. Although the orbital energies are influenced by the detailed distribution of 
charge in the cluster, they should reflect the effective charges on the atoms. They are a 
measure, through Koopman’s theorem, of the ionization potential, IP. The Is IP for a 
negatively charged atom is smaller than that for a neutral or positive one. The changes 
among the e ls ’s shown in Table II are consistent with the differences in the Mulliken 
population analysis charges. In particular, the central atom e ls in AIjq is larger (less 
negative) than the edge atom £ ls’ s - For MgjQ and AIjq, we also examined the 
wavefunctions for states very close in energy to those of Eq. (1). The behavior of the 
populations and the € ls ’s is similar to that shown in Table II. All our results lead to the 
central atom in Al 10 being significantly negatively charged while in Mg 10 or Be^g it is 
essentially neutral. 

We also examined 13 atom clusters designed to simulate the bulk environment for a 
central atom by surrounding it with its 12 nearest neighbors. For the hep systems. Be and 
Mg, the cluster was Xj 3(3,7,3) and we obtained SCF solutions for the lowest D 3 h closed 
shell state. For fee Al, the cluster was Xj 3 ( 4 , 5 , 4 ) and the SCF solution was for the lowest 
doublet state with a single open shell in O h symmetry. As for X 10 , the z axis is defined 
normal to the cluster layer. The valence level populations of the central atom for these X 13 
systems is given in Table III where the populations are divided into s, P x +Py and p z 
character. The symmetry of the clusters is such that p x =p y . The X 13 populations are 
compared to those for the X| 0 clusters. As in the case of the two-layer cluster, the central 
Al has a far larger charge than do the Be or Mg central atoms; in Al 13 (3,7,3), it is essentially 
Al~. For Be and Mg, increasing the number of neighbors also increases the charge on the 


central atom, but by much less than in the case of Al. 



The effect of the addition to the cluster of the atoms which are second nearest 


neighbors of the central atom is also of interest. However, this would lead to rather large 
clusters if we extend the X 10 or Xj 3 systems. Therefore, we chose to investigate this effect 
with clusters containing only one layer of atoms. For Be, we considered only Be 7 containing 
the nearest neighbors of a central atom on the (1000) face. For Al, we considered the (100) 
face and used Alj and A1 j 3 clusters. The Alj cluster contains a central atom and its four 
nearest neighbors; for Al 13 , we added the 8 second nearest neighbors of the central atom to 
Alj. We obtained SCF solutions for the lowest closed shell configuration for Be 7 and for the 
lowest single open shell configurations for Al 5 and A1j 3 . The central atom populations for 
these clusters are also given in Table III. 

Compared to the X 13 populations, the charge on the central atom in Al s is reduced by 
0.21 electrons; in Be 7 , it is reduced by 0.11 electrons. The central Al atom in Al 5 is still 
significantly negative while all the atoms, central and edge, in Be 7 are nearly neutral. The 
most interesting result comes when we consider the central atom charge in Al 13 which 
contains both first and second nearest neighbors. The change is reduced by 0.69 electrons 
from that in Al 5 ; the central atom is now nearly neutral! The central atom in Al 13 has 
virtually no overlap population with the modestly distant second nearest neighbors. 

Moreover, the overlap population between the central atom and its nearest neighbors is 
almost the same for Al 5 and Al 13 . Hence, clearly, the electron distribution about the central 
atom has changed dramatically when its second nearest neighbors are added to the cluster. 
The change in the population analysis charge from -0.78e in Al 5 to -0.09e in Al 13 is not a 
consequence of artifacts in the analysis. There is also a large change in the charge 
distributioirfor the central atom’s nearest neighbors between Al 5 and Al 13 . They were 
positively charged 0.20 in the small cluster and are negatively charged by 0.37 electrons in 



A 1 j 3 - The entire view of how bonding would take place could be altered by changes of this 
magnitude. 

As for Al 10 (7,3), we found two low-lying states for Al 13 (13,0). For Al 13 , unlike the 
results for AIjq, the charge distributions for the two states are fairly different. The central 
A1 atom is more nearly neutral for both Al 13 states than it was for AI 5 . However, the 
second state, higher in energy than that discussed above by only 0.04 eV, has a charge on 
the central A1 atom of 13.31. Thus, the addition of second nearest neighbors is not sufficient 
to insure that all the low lying states have a near neutral central atom. We conclude that in 
addition to the total energy of various states, the charge on the metal atoms directly involved 
in forming the chemisorption bond should also be used as a criterion for choosing the state to 
be used in studying adsorption. For the simulation of a metal surface, these atoms should be 
nearly neutral. 

Another interesting question concerns the extent of the hybridization of the central, 
adsorption site, atoms in the clusters. It is known, for example, that in many A1 compounds, 
a 3s electron is promoted to 3p leading to trivalent Al. We find that the central atom in the 
three-dimensional Al 13 (4,5,4) cluster has a 3s population of 1.12e (see Table III) consistent 
with a trivalent atom. We find similar central atom hybridization for the two layer A1 jq( 7,3) 
and both the one layer clusters, Alj^.O) and Al 13 (13,0). It would appear that the first 
nearest neighbors are sufficient to lead to a substantial 3s to 3p promotion. Previously , 3 we 
observed that there was a very large change in the Mg hybridization, as measured by the 
number of 3p electrons, between tetrahedral Mg 4 , 0.22, and Mg 13 (3,7,3), 0.85. We note 
that in Mgjg(7,3) the total 3p population is 0.54e; less than Mg 13 (3,7,3) but considerably 
larger than tetrahedral Mg 4 . The ns to np promotion for the central atom in Al 10 and Bej 0 



is also less than for A1 j 3(4,5,4) or Bei3(3,7,3), respectively. However, in all three cases 
there is substantial s to p promotion for the central atom of the X 10 clusters. 

Ideally one would like to use a cluster sufficiently large so that all problems of charging 
and hybridization arising from its finite size are eliminated; this is not always possible and 
one must find the correct balance between computer time, finite cluster effects and the 
accuracy needed to resolve the problem. In the past we have found is quite useful to use 
small clusters for the preliminary calculations, thus reducing the number of expensive large 
cluster calculations. The question naturally arises whether the results obtained with small 
clusters are meaningful. The answer depends, at least in part, on the strength of the 
interaction and the quality of the information desired. For weak interactions, charge buildup 
of the magnitude found for Al 10 may lead to seriously questionable results; for very strong 
interactions, it might add some uncertainty but useful information might still be gained. It 
would also appear reasonable to investigate initially a large bare cluster. Then to proceed 
with smaller clusters which have a similar charge distribution as the large cluster for the 
atoms which will be directly involved in bonding to an adsorbate. Preliminary results appear 
to support this view. Also, our recent work shows that sometimes the charge distribution for 
excited states of the smaller clusters are closer to that of the large cluster than is that for the 
small cluster ground state. Thus, it may be appropriate to use these excited states for the 
study of chemisorption. In this way, it should very often be possible to find a small cluster 
with reasonable charge distributions. However, it should be cautioned that there may be 
problems associated with performing SCF calculations on excited states if the addition of the 
adsorbate lowers the point group symmetry of the cluster. 

Our previous work for the chemisorption of H/Be(1000) provides guides for the size of 
a cluster needed to reduce the artifacts, or errors, associated with using a finite cluster to 



model a semi-infinite surface. However, it is clear that the magnitude of finite cluster effects 
varies from material to material and is even different for different faces of a crystal. 
Therefore, the question of the appropriate cluster size must be investigated separately for 


each chemisorption problem. 
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Table I. Summary of the height of H above the surface, r e in A, vibrational frequency for 
motion normal to the surface, « e in cm -1 , and adsorption site Be atom populations for 
H/Be(1000) at the open site. The Be populations are given for one of the three equivalent 
• adsorption site atoms. They are divided into 2s, 2p x , 2p y , 2p z and 2p=2p x +2p y +2p z 

character. The cluster surface is taken as the xy plane. 


Cluster 

2s 

2p x +2p y 

2 Pz 

2p 

r e 

"e 

Be 3 (3,0) 

19 

0.45 



1.21 

1230 

Be 6 (6,0) 

wtm 

0.86 

0.38 

1.24 

1.12 

1150 

Be 6 (3,3) 

wSM 

0.83 

0.41 

1.25 

1.05 

1200 

Be 9 (6,3) 

0.77 

0.78 

0.52 

1.30 

1.05 

1160 

Be 22 (14,8) 

0.62 

— 

— 

1.41 

0.94 

1230 

Be 36 (14,8,14) 

0.54 

1.05 

0.42 

1.48 

0.90 

i 

1200 
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Table n. Atomic charges and Is orbital energies, e ls in hartrees, for the Xjq( 7,3) clusters 
for Be, Mg and AL Results are given for the three distinct types of atoms; central, first layer 
edge and second layer. 


Atom 

Type 

Be 

Charge e ls 

Mg 

Charge e ls 

A1 

Charge e ls 

Central 

3.93 

-4.6203 

11.91 

-49.0656 

13.63 

-58.1445 

Edge 

Second 

3.98 

-4.6108 

12.00 

-49.0695 

12.90 

-58.2005 

Layer 

4.07 

-4.5847 

12.02 

-49.0366 

12.99 

-58.2085 
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Table HI. Central atom populations for one layer and two layer surface clusters and for 
three layer bulk clusters for Be, Mg and Al. The total population is divided into s, p x , p y 
and p z character. The xy plane is the taken as the plane of the layers of the clusters. 




Three Layer 

Two Layer 

One Layer Clusters 



*13 

X 10 

First Neighbors 

Second Neighbors 


'"total 

4.09 

3.93 

3.98 


1 

i 2s 

0.65 

0.79 

0.42 


Be< 

2Pz 

0.82 

0.34 

0.86 

— 


l 2 P X + 2Py 

0.62 

0.80 

0.70 

— 

I 

/total 

12.11 

11.91 

— 

— 

i 

Mg< 

|3s 

1.23 

1.36 

— 

— 

PPz 

0.33 

0.12 

— 

— 


(3P X +3 Py 

0.50 

0.42 

— — 

— 

| 

'"total 

13.99 

13.63 

13.78 

13.09 


1 3s 

1.12 

1.21 

1.29 

1.29 

All 

3P Z 

0.92 

0.70 

0.62 

0.37 


^3p x + 3p y 

1.84 

1.70 

1.88 

1.42 
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ABSTRACT: Self-consistent-field and correlated wave functions have been computed for the 
ground and for several low-lying states of CuO. The ground state is X 2 II and the lowest 
excited state, at ~8,000 cm -1 above X 2 II, is a previously unidentified 2 2 + state. The 
separation of these states is compared to that for the similar states of KO and is analysed in 
terms of integrals between orbitals of the separated free ions. We consider a classification of 
the states of the molecule based on states of Cu + and O' which leads to a division into 
manifolds of states arising from Cu + 3d*® and Cu + 3d^4s*. We predict that the states of 
the 3d 9 4s* manifold are 10,000-30,000 cm' 1 above -the ground state and assign the observed 
A 2 2 + state at 16,500 cm' 1 to this manifold. 


♦Present address: Polyatomics Research Institute, 1101 San Antonio Road, Suite 420, 
Mountain View,. California 94043. 

••Present address: NASA Ames Research Center, Moffett Field, California 94035. 



INTRODUCTION 


As a preliminary step for the study of the chemisorption of oxygen on copper, we have 
investigated the bonding in the diatomic molecule CuO for both the ground and the low lying 
excited states. We find that the two lowest states are ionic and arise from Cu + in a 3d 10 
configuration; these are the X 2 II ground state and a previously unidentified A 2 2 + state. We 
have also investigated additional ionic states which arise from Cu + in a 3d^4s* configuration. 
We conclude that the state which had been identified 1 as the first excited state at 
T e = 16,491 cm* 1 is likely to arise in this way. The separation between the X 2 II and A 2 2 + 
states is much larger for CuO, ~ 8,000 cm -1 , than for the similar ionic molecule KO 1,2 where 
it is only ~350 cm* 1 . This difference is analysed in terms of the overlap and nuclear 
attraction integrals between the anion and cation atomic orbitals. 

We consider a classification of states for CuO based on the properties of the separated 
atoms. The Cu atom has a 2 S(3d 1(1 4s 1 ) ground state; the first excited state, 3 3 D(3d^4s 1 ), is 
at 12,019 cm' 1 . For Cu + , the lowest state is ^Od 10 ) and the first excited state is 
3 D(3d 9 4s 1 ) at 22,648 cm' 1 . One expects CuO to be a very ionic molecule with essentially 
Cu + and O* character. The O* anion can, in C -v symmetry, be formed in two ways: aV 4 
(a hole) or a 2 ir 3 ( ir hole). When O* is combined with Cu + , ^Od 1 ®), a 2 II and a 2 2 + state 
can be formed. This is similar to the behavior found for the alkali monoxides, 2,4 LiO to 
CsO. We refer to these low-lying states as forming a manifold arising from 3d 10 . 

If O' is combined with the Cu + in the 3d 9 4s 1 excited configuration a large number of 
states can arise. These states can be either ionic or mixed (ionic and covalent) in character. 
For the first type, those with essentially only ionic character, three II, two 2 + , two A, one 
2' and one state come from the interaction of O" having either a a or a v hole in the 2p 
shell with Cu + having either a a, it or 8 hole in the 3d shell. The spin coupling of these 



states is doublet or quartet. We refer to this class of states as purely ionic states. We should 
note that the Mulliken population for either the SCF or the multi-configuration SCF 
(MCSCF) wavefunctions for these states shows a charge of ~0.7 and an overlap charge of 
~0.3. Thus, they are not purely ionic; however, we refer to them in this way to distinguish 
them from states which also have a covalent bond. This second type of states arises when a 
covalent bond is formed between the O' 2pa and Cu + 4s. This can lead to three states, 
2 2 + , 2 n, or 2 A, depending on the location of the d hole in Cu + (3d^4s 1 ). We call these 
mixed states since they have both a covalent a bond and an ionic contribution coming from 
Cu + O'. We refer to the solely ionic and mixed states as forming a manifold arising from 
3d^4s 1 (for short, the 3d 9 4s 1 manifold.) Given the separation between the 3d 1 ® and 3d 9 4s 1 
states of Cu + and the 3d 10 4s 1 and 3d 9 4s 2 states of Cu, it is reasonable to expect the Sd^s 1 
manifold of CuO to lie between 10,000 to 20,000 cm' 1 above the ground state. It is useful 
to consider this manifold as being formed from three groups of states. The first group 
comprises the three mixed states. The second group is composed of the ionic states for 
which O' has the configuration 2pcr 2 2pir 3 and the third group of states where O' is 
2p<x 1 2pir 4 . We will present results for representative states of CuO in each of these groups. 

For the CuO states discussed above, we have obtained LCAO SCF, MCSCF and single 
and double excitation configuration interaction, SDCI, wavefunctions. The basis set integral 
and SCF calculations were carried out using the MOLECULE-ALCHEMY program system. 3 
The MCSCF calculations were performed using the CASSCF methods and programs 
developed by Roos et alf* and the SDCI calculations using the graphical unitary group direct 
Cl program developed and implemented by Siegbahn. 2 

The organization of the paper is as follows. The details of the Gaussian basis sets are 
presented in Section II. In Sections III and IV, we discuss the nature and separation of the 


lowest, 3d 10 manifold, 2 II and 2 2 + states. We analyse and explain the role of the d 
electrons for the separation based on the overlap and nuclear attraction one-electron 
integrals. In Section V, we consider the nature of the excited CuO states that arise from the 
3d 9 4s l manifold. 

n. BASIS SET DETAILS 

Wavef unctions for CuO and KO have been computed in a basis of contracted Gaussian 
type orbitals. For both Cu and K, we have used the exponents optimized by Wachters 8 
including two functions which represent the 4p orbital. (Those for 3d i0 4s 1 , 2 S were used for 
Cu.) We have added a diffuse d function to this set. For Cu, we used the value determined 
by Hay^ and for K, we chose the exponent a=0.4. These basis functions were contracted 
slightly differently than recommended by Wachters 8 in order to have the largest flexibility for 
either ionic (cationic) or neutral character for the metal 3p orbital. The oxygen basis set is 
based on Dunning’s 10 double zeta contraction of Huzinaga’s primitive set to which we have 
added a diffuse p, a=0.059, to represent the O" character and a d, a=0.90. These added 
exponents have values close to those given by Dunning and Hay. 11 The final basis sets are: 
Cu(14sl Ip6d/8s6p4d), K(14sl lpld/8s6pld), and 0(9s6pld/4s3pld). 

III. STATES ARISING FROM 3d 10 MANIFOLD 

Symmetry and equivalence restricted SCF calculations were performed for the X 2 II and 
the lowest 2 2 + states. Cl wavef unctions consisting of all single and double excitations 
(SDCI) from the SCF reference were also obtained for these states. The SDCI excitation 
level is defined for a C xV reference state. Thus for the 2pir 3 X 2 II state, configurations 
which involve certain excitations of three spin orbitals are also included. In these SDCI 
calculations, excitations were permitted for only 17 electrons. The lowest 6 a and 2v SCF 
orbitals were frozen. These orbitals correspond to the Ar core of Cu and to the Ols. 


Excitations were not allowed into the highest five a and two v virtual SCF orbitals. For our 
basis, these orbitals are clearly functions which would correlate the frozen cores. 

The results of these calculations are summarized in Table I. The separation between 
the X 2 II and the A 2 2 + states is only slightly changed by the inclusion of correlation. The 
change in the separation between the SCF and Cl calculations is only 0.14 eV while the total 
correlation energy is ~10.8 eV. The small change is consistent with the view of an ionic 
system with small interatomic (molecular) correlation. The computed separation of only 
~8,000 cm' 1 is very different from the reported value 1 of 16,491 cm' 1 . This leads us to 
conclude that the state previously identified as the A 2 2 + state is not the first excited state. 

It is likely that it is a higher state arising from the 3d^4s 1 manifold, as we will discuss below. 

The computed equilibrium bond distance, r e , for the X 2 II state at the SCF level is 

O 

much longer than the experimental value of 1.72A. When correlation is included, the bond 
shortens. This is similar to other calculations on transition metal compounds, 12 ' 14 where 
correlation also shortens the bond. Even at the SDCI level, the agreement with experiment is 
poorer than one might have expected. However, the well is very shallow (see the vibrational 
frequencies given in Table I), and thus small errors due to the basis set or level of correlation 
can result in large errors in r e . Based on our calculated results, the experimentally 
unidentified 2 2 + state will have an r ft which is the same or slightly shorter than that for the 
X 2 n state and a slightly lower vibrational frequency, co e . 

We further investigated the correlation effects on the bond length by performing a 
multi-reference SDCI calculation for the 2 2 + state using SCF orbitals. Once again we 
correlated 17 electrons and used the same orbital space as in the single reference calculation, 
only now in addition to the SCF reference we also included all singles and doubles from a 


second reference configuration formed by a single excitation from the Cu3d<r to 02p<j. This 
second configuration is the only additional configuration to have a coefficient greater than 
0.05 in the single reference Cl; it has ionicity Cu +2 and O' 2 . The reference configurations 
were; 

(02s) 2 (02PCT) 1 (02pw) 4 (Cu3do) 2 (Cu3dir) 4 (Cu3d6) 4 - SCF 
(02s) 2 (02pa) 2 (02p*r) 4 (Cu3d a) 1 (Cu3dir) 4 (Cu3dfl) 4 - 3da to 2pa 

For other transition metal oxides CrO, MoO, NiO and PdO, we have seen that the inclusion 
of this type of configuration (nd to 2p excitations) leads to a considerable bond shortening. 
However for CuO it had a very small effect. The 2 2 + bond length shortened by only 
~0.02 a.u. upon going from single to multi-reference SDCI. 

We have also computed the dissociation energy, D e . It is known that for an ionic 
molecule 15 it is better to compute the D e relative to the ions, (Cu + and O') and to correct 
this value with the experimental atomic ionization potential and electron affinity. 16 We have 
computed the D e at both the SCF and Cl levels and find a D e (SCF) of 1.40 eV and a D e 
(Cl) of 1.88 eV. For a purely ionic system, the D e (SCF) and D e (Cl) are extremely 
similar, but in CuO they differ by 0.5 eV. This is consistent with the Mulliken populations 
which show some covalent character. However, the Cl value is much smaller that the 
experimentally cited value 1 of 2.79 eV. Even considering the errors which may arise due to 
the choice of basis set and level of the correlation treatment, it appears 2.79 eV is quite likely 
to be too large. We estimate that a better value for D e is in the range 2. 0-2. 3 eV. 

IV. ORIGIN OF THE 2 II- 2 Z + SEPARATION 

One important feature of CuO is the size of the lowest X 2 II-A 2 2 + separation 
(~7,200 cm' 1 ). The Cu + ion 1 S(3d 10 ) is a closed shell state as is K + ^Op 6 ); but in KQ 


the 2 n- 2 2+ separation 1,2 is only ~350 cm -1 . In order to analyse the origin of this different 
behavior, we calculated some molecular properties for KO and CuO using SCF orbitals for 
K + , Cu + and O'. (For O', the atomic wavefunctions were computed in C ooV symmetry. We 
used the orbitals from the 2pa 2 2pir 2 state.) The molecular properties were computed for a 
bond length near equilibrium; r=4.25 bohr for KO and r=3.27 bohr for CuO. The energies 
for both 2 I1 and 2 2 + states were evaluated for the atomic orbitals with the O" orbitals 
Schmidt orthogonalized to those of the metal ion (K + or Cu + ). We refer to the difference 
between this frozen orbital energy and the SCF energy as the relaxation energy. For KO, 
this relaxation energy is very small, ~5,000 cm' 1 . This is consistent with a highly ionic 
behavior which is characterized by a Mulliken gross charge of ~0.9e on K and a very small 
K-O overlap population of ~0.07. For CuO, the relaxation energies are somewhat larger, 

~ 13,000 cm* 1 . This is consistent with a less ionic gross charge of ~0.7 and a larger overlap 
population of ~0.3. The differential relaxation energies between the 2 II and 2 S + (given in 
Table II) are very similar for KO and for CuO; for both molecules it is rather small 
~ 1,800 cm' 1 . These calculations support the view that 2 2 + and 2 II states for both systems 
are dominately ionic. 

The origin of the separation of the 2 2 + and 2 II states can thus be analysed in terms of 
integrals between the orbitals of the separated ions. We have evaluated the overlap of the 
02pa and 02pir with the metal, Cu + or K + , orbitals and the one electron, kinetic energy 
plus nuclear attraction, integrals for the oxygen p orbitals denoted h(2p<r,2p<y) and 
h(2pir,2pw). The overlap integrals are a measure of the Pauli repulsion and the one electron 
integrals a measure of the nuclear attraction. To account for the shielding by the core 
electrons, the core potential was evaluated for the configurations Culs 2 ...3d 10 , Kls 2 ...3p 6 , 
and 01s 2 2s 2 and added to the one electron integrals. This potential is SNjUj-fcKj); 
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i=ls,2s,..., where Nj is the orbital occupation number and Jj and Kj are the Coulomb and 
exchange operators respectively. In effect, the metal atom is treated as a + 1 ion. These 
quantities are summarized in Table II. For KO, the difference between h(2pir,2prr) and 
h(2pa,2pa) is small, however the a-a overlap is 6 times the ir-ir overlap. Thus for KO, the 
repulsion dominates and the 2 2 + state which has one less a electron is the lower state. Also 
the 2 n state has a longer r e , 2 this lengthening reduces the overlap repulsion. The net 
behavior is that the separation of these states is very small. In the CuO case, the difference 
in the one electron integral is 4.4 times larger than in KO. The larger value for h(2pa,2pa) 
favors the 2 II state which has two 2pa electrons. The ratio of o-o to ir-v overlap is 2.5, 
down from 6 in KO. The lower a repulsion and the increase in the a nuclear attraction, 
h(2po,2pa), relative to ir leads to a 2 II ground state and a much larger 2 II- 2 2 + separation. 
Also since a replusion is not as dominant a feature for CuO, the r e for 2 II and 2 2 + are the 
same at the SCF level, while for KO, they differ by 0.1 2A. Similar effects have been used 
to explain the 2 2 + - 2 n separation for the alkali metal monoxides. 4 


# 


V. STATES ARISING FROM 3d 9 4s 2 MANIFOLD 

To investigate the nature of the excited state spectrum of CuO, and, in particular, to 
identify possible candidates for the experimentally observed 2 2 + state at -16,000 cm' 1 , we 
have performed SCF, MCSCF and SDCI calculations on several of the states in the 3d 9 4s 1 
manifold. Our interest is to obtain qualitative information on the nature and positioning of 
the states in this upper manifold. 

We expect the two ionic groups of this manifold. O' 2pa 1 2pir 4 and 2pa 2 2pir 2 , to be 
similar to the 3d 10 manifold, with the ionic O* 2pir hole states lower by -8,000 cm' 1 . Each 
group will, of course, be split into several states because of the three possible d holes (o, ir 
and 5) and the two different spin couplings (doublet and quartet). We have investigated 


various orientations of the d holes, and for both of the ionic groups of states, as well as for 
the mixed group of states, we find the splitting to be small. Furthermore, we would expect 
the splitting between these two ionic groups to be well computed at the SCF level as was the 
case for the 3d*® manifold. The third group, the mixed states, is different and it is not easy 
to predict their location relative to the two ionic groups. It can be considered as having one 
02pc electron, therefore, as an upper bound it would lie near the ionic 02pa hole state; 
however, it is lowered by the strength of the Cu(4s)-0(2p<r) covalent bond. The strength of 
this bond is not well computed at the SCF level. The minimum treatment required is a 
simple MCSCF accounting for the correlation of this a bond which is not present in the two 
ionic groups. Thus, a qualitatively accurate splitting among the groups of the 3d^4s* 
manifold should be given by SCF calculations for the two ionic groups and a simple MCSCF 
for the mixed group. 

In addition to the splitting within this manifold, we would like to position the 3d 9 4s* 
manifold with respect to the 3d 10 manifold. It is known that at the SCF level, the 
Cu^D(3d 9 4s 2 )- 2 S(3d*°4s*) and the Cu + *S(3d 1 ®)- 2 D(3d^4s 1 ) splittings are in poor 
agreement with experiment (as shown in Table III) and we expect similar errors to arise in 
the CuO molecule. Since it is possible to clearly identify the molecular states as arising from 
3d^ or 3d 10 , we simply shift the computed results by a correction factor to account for the 
intrinsic errors in the SCF description of Cu and Cu + . The correction (11,460 cm' 1 ) is 
computed to be the weighted average of the SCF error in Cu (~9,000 cm' 1 ) and Cu + 

(~ 12,000 cm' 1 ) based on a charge of 0.75 which is approximately that computed for CuO at 
the SCF level. Since the errors of Cu and Cu + are reasonably similar, this correction would 
not change greatly if a somewhat different choice of ionicity were used. 
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At the Cl level, the error in the Cu and Cu + splittings has been considerably reduced, 
see Table III; thus for CuO, we expect that the uncorrected Cl results for the vertical 
separations will be quite similar to the corrected SCF (or MCSCF) results. The Cl 
calculations serve both to refine our separations and as a check on the simple model on 
which the SCF/MCSCF results are based. 


* 


•> 


For the states of the 3d 9 4s 1 manifold, we carried out all calculations at 
R(Cu-0)=1.73A. the experimental r e of the X 2 II state. 1 The r e for all observed states is 
very similar 1 and hence we felt that the excited states could be computed at one point. 

The three different states of the mixed group were evaluated at the SCF level and the 
results are summarized in Table IV. We find the splitting among them is only 0.4 eV. We 
also found similar splittings, upon changing the orientation of the d hole, in the two ionic 
groups, thus for further calculations we considered only a single 2 II state for each group as 
representative of all of the states in the group. The 2 II symmetry for the representative state 
was chosen for computational convenience. 

The MCSCF wavefunctions were obtained using the CASSCF technique. 6 In this 
approach the orbitals are divided into inactive orbitals, which are doubly occupied in all 
configurations and active orbitals which are occupied with all possible symmetry allowed 
distributions of the active electrons. The same MCSCF was performed for all three 2 II 
states. It consists of 5 active electrons with 2 active a and 1 active it orbital. We have also 
imposed the restriction of 2 electrons in the a space and 3 in the v This leads to 4 
configurations. (Symmetry and equivalence restrictions were imposed for the it but not for 
the S orbitals. We tested the importance of this constraint with the SCF wavefunctions and 




found it to lead to changes in the energy of less than 0.2 millihartrees.) For the mixed state. 


the active a orbitals are the bonding and anti-bonding combinations of Cu, dominantly 4s 
and 02p a orbitals; the active it orbital is essentially Cu3dw. For the ionic 02p a hole state 
the active a orbitals are essentially Cu4s and 02p and the active it is 3d it. For the ionic 
02p it hole state the active orbitals are Cu4so, Cu3da and 02pir. This choice of active 
orbitals and electrons leads to a large correlation effect, ~ 14,000 cm' 1 , for the mixed state 
(see Table V) which is typical of a covalent bond. For the two ionic states we have 
effectively performed an SCF calculation since a single CSF (configuration state function) 
has a weight of >0.99. The corrected vertical excitation energies from these calculations, A ; 
in Table V are ~ 15, 000 cm' 1 for the mixed and ionic 02pir hole states, the 02p a hole state 
is much higher at ~27,000 cm' 1 . 

The SDCI calculations for these states were performed with the same number of 
deleted virtual and frozen, doubly occupied, core orbitals as for the calculations on the 3d 10 
manifold. The orbitals used were those from the MCSCF calculations described above. The 
excitations were made from only one. component of the degenerate 2 II state. For the two 
ionic states, the single and double excitations were made from a single configuration, these 
CSF’s are shown in Table V. The mixed state excitations were made from two CSF’s, one is 
shown in Table V and the second is formed by replacing the bonding (4s+02po) 2 with the 
anti-bonding (4s-02pa) 2 in this CSF. The Cl vertical excitation energies, A in Table V, are 
given as calculated without correction. 

The difference between the MCSCF and Cl results for the separation of the ionic 
states and for the Cl A and the corrected MCSCF A ; are rather small. This is consistent 
with the small difference between the SCF and Cl separations for the 3d 1 ® X 2 I1 and A 2 2 + 
states. It also follows because the Cl leads to a very small mixing between 3d 10 and 3d 9 4s 1 
configurations. There is some change between the MCSCF and Cl results for the mixed 
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state, the Cl A is ~3,000 cm' 1 less tham the MCSCF A ; . Although the MCSCF 
wavefunction is able to qualitatively describe the bonding pair, the extensive correlation of 
the SDCI treatment is needed to position this state more accurately. This Cl treatment 
brings the mixed and the ionic 02pir hole states virtually on top of each other. We conclude 
that the state, observed 1 at ~ 15,000 cm' 1 is either a mixed state or an ionic 02pir hole state. 
The ionic 02p a hole state is too high to be considered. The separation between the 02p it 
hole and 02pa hole ionic states is larger than in the 3d 10 manifold. We assume that this is 
because the 02pir hole state has only one 3 da electron thus reducing the a overlap and 
further enhancing the nuclear attraction preference for the 02pir hole over the 02pa hole 
state. 


e 


We should note that our MCSCF and Cl calculations for this manifold do not satisfy 
McDonald’s theorem. Both the mixed and ionic 02pir hole states were the lowest roots of 
the CASSCF calculation, while the ionic 02pa hole state was optimized as a second root. In 
all cases, the lowest root of the Cl was selected. The fact that the lowest roots are 3d 9 4s J 
like and not 3d 10 like follows simply because the orbitals for the appropriate states of the Cu 
atom and the Cu + ion are different. If the one configuration energy of Cu 2 D(3d 9 4s 2 ) is 
evaluated using the SCF orbitals determined for Cu 2 S(3d 10 4s 1 ), the separation between 
them is 4.70 eV. If separate SCF calculations are performed for the 2 S and 2 D states, the 
separation is 0.37 eV. A Cl calculation based on a common set of orbitals for these states 
would have to be very large to describe these large orbital relaxation effects. Our Cl 
calculations (~ 100,000 configuration) were not yet large enough to account for them. Thus 
each SDCI calculation was based on individually optimized orbitals from the appropriate 
MCSCF wavefunctions. This accounts for the relaxation simply and in a compact fashion. 
Since the Cl’s for the 3d 10 manifold do not appear to mix in much 3d 9 4s 1 character, we feel 



the fact that there is no lower 3d 1 ® like root in these calculations is not a severe problem. 
That the ionic 02pir hole and the mixed states are virtually identical in energy and that each 
is the lowest root is a potential problem. However, we feel that we have shown that the 
state observed at ~ 15, 000 cm' 1 is due to either a mixed or an ionic 02pir hole state and that 
the ionic 02pa hole state is much higher up. It is certainly true that a more extended 
calculation might lead to a mixing, and hence somewhat separate the mixed and 02pir hole 
states. However, this is outside our goal of qualitative understanding of the 3d^4s 1 manifold. 

VL CONCLUSIONS 

There is an unobserved 2 S + state at ~8000 cm' 1 arising out the the Cu + 3d 10 
manifold. The previously identified A 2 2 + is probably a state arising from the Cu + 3d 9 4s 1 
manifold. While we have performed calculations on 2 II states, we expect these to lie close to 
the 2 2 + states, since our tests on the Sd^s 1 manifold indicate that the energy of the state is 
reasonably independent of the position of the d hole. We found both the mixed state and the 
ionic state with an 02pw hole to be positioned at ~ 15,000 cm' 1 , and considered these as two 
candidates for the experimentally observed state. The third type of state in the 3d 9 4s 1 
manifold, Ospo hole, is at ~30,000 cm' 1 . 

A comparison with KO shows that while the d’s are not directly involved in the 
chemical bonding, their presence leads to a change in the ratio of the a one electron integrals 
to those of v symmetry. The a are closer to the it overlap integrals for CuO than for KO. 
The a one electron integral is also considerably smaller than the v integral. Thus, for CuO, 
both the nuclear attraction significantly favors the 2 II state and the Pauli repulsion bias 


against it is reduced. 
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TABLE I 

® The SCF and Cl results for the X 2 II and A 2 2 + states of the 3d 10 manifold. The properties 

given are the equilibrium Cu-O bond distance, r e , the vibrational frequency, « e , and the 
dissociation energy, D e . Also given are the total energy, E, the energetic separation of 
the states, A, and the correlation energy obtained for the SDCI wavefunctions. The energies 

are those computed at the appropriate r e ’s. 

m 


/ 


€> 




2 n 

2 2 + 

A (cm* 1 ) 

SCF 

r e (A) 

1.875 

1.875 



E(a.u.) 

-1713.59823 

-1713.56546 

7190 


o> e (cm -1 ) 

527. 

488. 



D e (eV) 

1.40 

— 


Cl 

r e 

1.82 

1.81 



E 

-1713.99679 

-1713.95889 

8320 


w e 

545. 

534. 



D e (eV) 

1.88 




Ecorr (eV) 

10.85 

10.71 


Experiment 3 

r e 

1.724 

— 



w e 

640.17 

— 



D e (eV) 

2.79 

— 



a See Ref. 1. 
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TABLE D 


Tabulation of the one electron properties for KO and CuO. For the Mulliken populations 
the gross charge is followed by the overlap charge in parentheses. The energies of the 2 II 
and ^2 + molecular states are also given; the SCF energy is compared to the frozen orbital 
energy obtained with the orbitals of the separated ions. The 2 II and 2 2 + separation, A, is 

computed in both ways. 



CuO r=3.27 a.u. 

KO r=4.25 a.u. 


Muliken Populations 


2 n O population 

0.719 (0.295) 

-0.909 (0.071) 

2 2 + 0 population 

0.771 (0.263) 

-0.923 (0.004) 


Overlap Integrals 


<02pa 1 Cu3da> 

0.091 


<02po | Cu3pu> 

0.081 <02p<r | K3p<r> 

0.149 

<02p a I Cu3s> 

0.126 <02pu|K3s> 

0.142 

<02pir | Cu3dir> 

0.088 


<02pir | Cu3pir> 

0.030 <02pir | K3pir> 

0.050 


One Electon Integrals (Hartrees) 


h(2pa,2pa) 

-3.0726 

-2.9108 

h(2pu\2pir) 

-2.9554 

-2.8842 


Energies 


SCF 2 n (Hartrees) 

-1713.59147 

- 673.94486 

Frozen orb. 2 II 

-1713.53085 

- 673.92844 

Relaxation (cm' 1 ) 

13303. 

3603. 

SCF 2 2 + (Hartrees) 

-1713.55972 

-673.94679 

Frozen orb. 2 2 + 

-1713.50738 

-673.92211 

Relaxation (cm' 1 ) 

11486. 

5416. 

ASCF (cm- 1 ) 

7089. 

425. 

AFrozen (cm' 1 ) 

5151. 

1390. 






* 



TABLE HI 


Summary of Cu and Cu + calculations. 

The separation. A, is in wavenumbers and the total energy in Hartrees. 



E(3d 9 4s 2 ) 

Cu 

E(3d 10 4s 1 ) 

A 

EOd^s 1 ) 

Cu + 

E(3d 10 ) 

A 

NHF a 

-1638.9505 

-1638.9662 

3450 

-1638.6814 

-1638.7287 

10380 

SCF 

-1638.7611 

-1638.7749 

3033 

-1638.4915 

-1638.5393 

10490 

SDCI 

-1638.9775 

-1638.0248 

10374 

-1638.6741 

-1638.7729 

21690 

EXP b 



12020 



22649 


a Numerical Hartree Fock. 
b Ref. 7. 
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TABLE IV 

SCF energies (in Hartrees) for the three mixed states arising from 
the 3d^4s* manifold, calculated at r=3.27 a.u. 


State 

Valence Configuration 

Energy 

2 Z + 

(4s+02pa) 2 3da I 3dir 4 3dfi 4 02pir 4 

- 1717.50912 

2 n 

(4s+02pa) 2 3d<7 2 3dir 3 3d5 4 02pir 4 

-1713.50494 

2 a 

(4s+02po) 2 3da 2 3d7r 4 3d6 3 02pir 4 

-1713.49360 




TABLE V 


SCF, MCSCF and SDCI energies (in Hartrees) and vertical excitation energies, A, (in cm* 1 ) 
for the three types of 2 II states arising from the 3d 9 4s 1 manifold and the 2 I1 state arising 
from the 3d 10 manifold. All states are calculated at r=3.27 a.u. For the Cl, the weight of 
the SCF configuration is given in parentheses. SCF and MCSCF vertical excitation energies 
corrected for the 3d 1 ®-3d 9 4s 1 separation error, A r , are also reported. 


Energy 


A 


A' 


A. 3d 10 State, X 2 ri-3da 2 3dir 4 3dS 4 02p<7 2 02pir 3 


SCF 

-1713.59823 



Cl 

-1713.99679 (0.91) 

0 

0 

B. 3d 9 4s 1 Mixed State-(4s+02pa) 2 3da 2 3dir 3 3d$ 4 02pir 4 



SCF 

-1713.50493 

20,480 

31,940 

MCSCF 

-1713.57022 (0.82) 

6,150 

17,610 

Cl 

-1713.9325 (0.85) 

14,110 


C. 3d 9 4s 1 02p<7 

Hole-4s 1 3da 2 3 dir 3 3 d5 4 02p<r 1 02pir 4 



MCSCF 

-1713.52615 (0.99) 

15,820 

27,280 

Cl 

-1713.8730 (0.91) 

27,170 


D. 3d 9 4s 1 02p7r 

Hole-4s 1 3d<r 1 3dir 4 3dfi 4 02pa 2 02pir 3 



MCSCF 

-1713.58735 (0.99) 

2,390 

13,850 

Cl 

-1713.9322 (0.91) 

14,180 
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Background and aim of the work 


The original MOLECULE integral program was written in 1971-1972 
/ 1 / . With minor modifications /2/, it has been used by quite a 
number of scientists around the world, and it has been implement- 
ed on many different computer systems. On the scalar machines 
available in the 70 's - and for the basis sets generally consid- 
ered adequate at that time - it performed rather well. The fact 
that it produces integrals over a symmetry-adapted basis still 
makes it an attractive alternative for many types of applica- 
tions . 


The ongoing development of vector-oriented computers and the 
increasing use of large basis sets with high Jl-values places new 
demands on an efficient integral program. In order to run effici- 
ently the code should be organized such that it benefits from the 
specific features of modern hardware, e.g. parallelism, pipe-lin- 
ing and special vector -oriented instruction sets. This is not the 
case with the old version of MOLECULE, and accordingly it doesn't 
run efficiently on supercomputers such as the CRAY-1 or CYBER- 
205. 


The aim of the present work has been to produce an integral code 
suited for modern, vector-oriented hardware, in particular for 
the CRAY series of computers. It should be capable of handling 
basis functions with high angular quantum numbers in an efficient 
way. The code is assumed to replace the old MOLECULE program, and 
a reasonably compatible in- and output has therefore been a 


desire. 


2 . 


General considerations 


The term " vectorization" is normally used for the process of 
adapting computer code to vector-oriented hardware. This usually 
amounts to modifying the code to make it benefit from the special 
vector-handling facilities available on modern computers. In a 
vectorized program, a major portion of the computational work 
should be performed in simple loops of sufficient length to 
utilize these facilities in an optimum way. 

In many cases, the vector structure is inherent in the algorithms 
and the original code. Vectorization then amounts to cleaning up 
the innermost loops, removing IF-tests and subroutine calls, 
minimizing I/O, etc. 

In other cases, however, the vector structure is not evident from 
the algorithm. The kernel of the code may be scalar, highly 
branched or contain loops of insufficient length. The execution 
is repeated through long outer loops over large portions of code, 
which do not invoke the use of special vector instructions. 

The two situations briefly described above are often referred to 
as "intrinsic" and "extrinsic" vectorization cases. Clearly the 
extrinsic vectorization is by far the more difficult one, as it 
requires changing the structure of the entire code - sometimes 
also of the algorithms. 

Many tasks involving manipulation of vectors and matrices can be 
performed in alternative ways. If the algorithm contains nested 



loops there is usually a freedom of choice as to the order of 
these loops. The optimum choice may be dependent on input para- 
meters, e.g. on vector lengths, but also on the type of hardware 
used. The task of multiplying two rectangular matrices may serve 
as an example. In the triple-nested loop constituting the kernel 
of such an operation, one would normally choose the longest loop 
as the innermost one in order to obtain efficient vector perform- 
ance. In some cases, however, that would amount to evaluating 
each element in the result matrix as a dot product. That scheme, 
while allowing peak performance on an FPS164, would not vectorize 
at all on true vector machines like the CRAY or CYBER-205. Of the 
two other possibilities, only one would perform well on a 
CYBER-205 due to problems with a vector stride larger than one. 
The comparison is further complicated by the fact that the latter 
two schemes may exploit the sparseness of the initial matrices, 
which the scalar product does not. 

Evidently, it is not always possible to design a straightforward 
algorithm which is both optimized and portable. Our approach to 
the problem has been to code routines for such standard elements 
of computation, which involve several alternative possibilities. 
The final choice, made at execution time, is dependent on the 
loop lengths, and, in addition, on parameters in the code which 
may be optimized at each installation, depending on the hardware 
used and the type of problem for which optimum performance is 


desired . 


3 . Comparison with the old code 


The original MOLECULE code is a typical case for extrinsic 
vectorization. The innermost, easily recognizable loop structure 
is the nest of four loops over primitive GTOs which together make 
up one integral over contracted functions. Inside these loops 
there are some 1000-1500 lines of code. Extensive testing -and 
branching is performed to take advantage of the simplifications 
due to low quantum numbers, planar molecular fragments, etc. 

Calls to external routines also hinder the vectorization. 

The integrals over GTOs are evaluated in batches, made up by all 
integrals within a quadruple of GTSs (see Appendix A for 
definitions). This saves recomputation of several quantities, but 
the loops generating these batches are not well suited for 
vectorization. To get loops of sufficient length it was found 
necessary to change the entire structure of the program. 

In the new code, all GTOs on a center with a given 1- value are 
considered as a group ( GGTO ) . Certain compromises and 
modifications of that approach may be necessary due to finite 
memory size, as discussed further below. The innermost loops in 
the integral code is made over groups involving all integrals 
over functions in four such GGTOs. The vector length for those 
operations will become = n^»n 2 *n 2 *n 4 ? n^ being the number of 
GTOs in the i:th group. In a very schematic way, one may view the 
new procedure as derived from the old one by replacing every 
scalar operation in the kernel by a vector loop of length N^. 


4 . Memory management 


The dimensioning of arrays in an integral code requires some 
special attention. Different arrays may become very large in 
different situations, depending on the quantum numbers, degree of 
symmetry, number of primitive functions in one contracted, etc. 
Dimensioning each array to the largest size attained in routine 
applications would quickly blow up the memory of most computers. 

On the other hand, it is quite unlikely that these arrays would 
all be large at the same instant of execution. (A molecule with 
eight symmetry-related atoms, each with five or six f-functions 
contracted to one is not a likely candidate for routine calcula- 
tions). Our approach has therefore been to use a truly dynamic 
allocation of meory. All large arrays are stored in two main 
working areas, and have variable dimensions which are passed in 
the subroutine calls. The array boundaries are allowed to change 
at execution time. With this arrangement each array takes only 
the space needed at any particular instant. In addition, tempor- 
ary arrays used at different stages of the calculation are allow- 
ed to overlap. 

However, even with maximum packing of arrays at any given 
instant, memory limitations may be encountered. This is rather 
easily handled for the case of segmented contraction. In that 
case, the loops over primitive integrals may be divided into 
several passes, using a "strip-mining" technique. In each pass, 
one computes only the portion of primitive integrals which can be 
easily handled in the available memory. It is trivial to sum 



these integrals as they arrive, in order to make up a contracted 
integral, and this step is hardly affected at all by the strip- 
mining. 

For the general contraction case, however, one would need all the 
primitive integrals before starting the transformation, but 
clearly this is not always feasible. The situation is further 
complicated as the program takes certain shortcuts whenever two 
or more of the GGTOs involved are the same. Apparently, the per- 
mutational symmetry present in the full integral set for such 
cases will be obscured if only part of the set is considered at a 
time. 

We have introduced a rather reasonable compromise in order to 
tackle this problem. In the segmented contraction case the effec- 
tive vector length is n^»n2*M; M ranging between 1 and n^*n^ 
depending on the memory available. For general contraction, how- 
ever, the choice of M had to be restricted to the values 1, n3 or 
n^*n^. (The n^:s are defined in Section 3). 


5. 


Computational simplifications 


In any computational work, it is customary before time-consuming 
sub-tasks are performed, to test whether the results of the sub- 
task will significantly affect the final result. 

When evaluating molecular integrals there is an almost unlimited 
number of special cases which may be treated with fewer numerical 
operations than the general case, and which, accordingly, are 
candidates for such tests. Simplifications may arise due to 
n-center cases with n<4, for integrals involving linear or planar 
fragments, for low angular quantum numbers, or because small 
multiplicative prefactors make the evaluation unnecessary. Sim- 
plifications also arise due to global - or local - symmetry (of 
which linear and planar fragments are examples ) . 

The extent to which pre-testing is worthwhile depends on the 
efficiency with which the test can be performed, the hit-rate of 
such a test, the complexity of the operations to be avoided and 
on possible degradation of program efficiency as a result of such 
testing. 


In the old program, integrals are computed in batches associated 
with a quadruplet of shells. Within a batch, the product of 
radial overlaps 


T , = S S, 

|ivX.o |iv \a 


(1) 


is a common factor which largely determines the magnitude of all 


integrals in the batch. By comparing this factor against a 
threshold. 


a single test may eliminate the evaluation of an entire block. 
With that scheme used in the new code, the testing would have 
occurred in the innermost loops over GTOs, and it would have 
inhibited vector ization completely. A different approach has 
therefore been taken: 

Initially, the sub-ranges of elements T > e are determined 

Hvya— 

for all T within the group to be processed. T, and all other 

4-index quantities necessary for the integral evaluation are then 
compressed into smaller arrays, based on those sub-ranges. The 
operation is nearly identical to the Q8VCMPRS on the CYBER-205, 
but since few other vendors offer that possibility we have coded 
it in Fortran. After compressing all arrays the integrals are 
evaluated, usually a smaller number than the full set. In 
cases of general contraction it is necessary to return to the 
full, ordered set of integrals before transforming to the 
contracted basis, This is achieved by expanding the array with 
the inverse operation, which again is very similar to the Q8VXPND 
on the CYBER-205. 

The old program made extensive use of special code for low quan- 
tum numbers. This led to substantial savings for certain types of 
integrals, as compared to the general procedure. With a vector- 
ized code the gain may be less, as the overhead in the general 


scheme is amortized over a large number of integrals. Further- 
more, in most time-critical applications today integrals with low 
quantum numbers play a minor role, anyway. We have therefore 
chosen to develop a first version without special routines, apart 
from some specific coding for ssss-type integrals. However, a 
vectorized code for low quantum numbers is in preparation, and 
will be implemented in the future. 

The simplifications due to local - or global - symmetry are 
treated outside the innermost loops, and may therefore be 
retained in the new code without degradation of performance. The 
situation when two or more GGTOs are the same is handled by fold- 
ing the arrays of radial overlap ” s ^ v and/or in ( 1 ) - when- 

ever they are symmetric. The resulting zeros in T ^ will then 
be removed by the COMPRESS feature of the program, which in turn 
enhances the speed. 



6 . 


Contraction 


The segmented contraction scheme (i.e. any primitive function 
contributing to only one contracted function) se4ems to present a 
problem on modern vector hardware, in the sense that it is 
diffricult to obtain efficient vectorization. The general 
contraction scheme / 3,4/ is more favourable in this respect; 
however, care must still be taken in order to achieve optimum 
performance. Apparently, a four-index n 5 transformation scheme 
involves the minimum number of operations. On the other hand, the 
index range in each step of the procedure will be too short to 
permit efficient vectorization. The transformation has therefore 
been arranged so that in any quarter-transformation the loops 
over the indices involved in the transformation are done outside 
the loops over the three passive indices. This has been achieved 
by performing the transformation steps 

[ap y 6 ] ♦ [|ip y 6 ] -*• [ m- v y6] -*■ [nv X 6 ] [nv \a] (3) 

using as index quadruples 

(a,P,y,6) , (p,y,6,u) » (y,6,|i,v) , 

( 6, n, v, \) , ( n, v, X, a) (4) 

respectively, where a,p,y,6 represent primitive functions, 
li,v,y,a contracted functions, and the index (a,p,y,6) denotes 
storage with the first index running fastest, etc. With n^ , n 2 » 
etc. primitives involved, and m^, m 2 , etc. contracted functions, 
a compressed passive index is defined from the first three 


indices of a quadruple in ( 4 ), say for (P,y,6,n) 


as 


p = Y + (6-1) n 3 + C^-l) n 3 n 4 


( 5 ) 


Given the running order of the various indices mentioned above, 
the location for the integrals [ p. p y 6 ] and [4V -y 6 ] involved in 
this particular quarter-transformation is given by 


P + (p-1 ) »n 2 ( 6 ) 

and 

p + ( v- 1 ) »n 3 *n 4 'mj^ ( 7 ) 


These integrals are accessed in a regular way, with the innermost 
loops over 4, y and 6 collapsed into a single loop over p. This 
general contraction scheme vectorizes efficiently on a computer 
where contiguous vectors are not necessary, such as the CRAY. The 
vector lengths range between m^*m2*m 3 and nj’n^n^. 


7. 


Specific difficulties encountered 


When the project was initiated, we anticipated to have access to 
the CRAY-1 in Sweden from late 1983 for program development and 
testing. The CRAY was delayed, however, and when it arrived we 
found the use of it extremely cumbersome due to bureaucratic 
obstacles of various kind. 

The code has therefore been largely developed on conventional 
mini- and mainframe computers, with or without an attached 
FPS164. Only occasionally we have taken it for testing on the 
CRAY* In retrospect, however, this turned out not to be as seri- 
ous a drawback as we had feared. Far more time was spent on a 
complication of an entirely different nature: 

The original plans were to produce a code that would offer the 
same possibilities as the old MOLECULE program, apart from a 
substantial speed-up due to vectorization over a large number of 
primitive integrals. However, with all these primitive integrals 
available simultaneously, the possibility of a general contrac- 
tion scheme" seemed rather obvious. This option of performing a 
general contraction was originally viewed as an extra little 
benefit from the new organization, to be exploited whenever poss- 
ible considering other restrictions. Only later did we realize 
the full power of the general contraction scheme. As it turns 
out, it provides quite a high-quality description of the atomic 
core even with a minimal basis representation. This seemed extre- 
mely attractive for large-scale post-Hartree Fock calculations on 
systems with heavy atoms, where it is important to keep the 


number of basis functions as low as possible. We therefore 
decided to make the general contraction scheme an integrated, 
fully supported feature of the new program. This brought some 
rather serious problems to the organization and design of the new 
code. Solving these, as discussed in Sections 4-6, without degra- 
ding the general performance of the program, required a rather 
delicate analysis and some elaborate programming, and delayed the 
project a lot. We feel, however, that the resulting code was 
worth the effort spent on it. In its present status it is capable 
of handling virtually any reasonable basis set within a memory of 
0. 5-1.0 MW, and most cases with a lot less than that. 


8 . 


Guide for Users 


This section is not intended to be the final program documenta- 
tion, to be used as a program manual (such a documentation is 
being produced and will be distributed soon). Here we just report 
some of the experience gained so far, along with some remarks, 
hints and warnings which may be helpful when trying to install 
and use the code in its present shape. 

As mentioned in Section 3, the primitive integrals are calculated 
in groups. Such a group may consist of integrals to be combined 
in a contraction scheme, either of segmented of general type. 
However, in most standard basis sets a large fraction of the 
basis functions consist of one primitive only. This would produce 
too short vectors to be of practical use if the above strategy 
were rigorously followed. Such "solitary" functions should 
therefore be grouped together and formally contracted by a 
one-to-one transformation. The program will recognize such a 
dummy contraction and take shortcuts to avoid unnecessary 
overhead. For optimum performance, the degree of such dummy 
contraction (in other words, whether e.g. six solitary functions 
should be combined to 

one, two or three groups) must be tuned to the hardware used. 

The size of the working areas are set in the main program. For 
best performance in large-scale applications, these should be set 
as large as conveniently allowed, considering other limitations 
and restrictions in the computing environment. Note, that for 
large applications it may be necessary to upgrade some of the 
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# 


small arrays. These are not dynamically allocated, and must be 
changed in a consistent way in all COMMON statements throughout 
the program. 

If the memory allocated should be insufficient, despite the fea- 
tures described in Section 4, the program will normally stop 
after a short message is issued. This procedure is not 100 % 
fool-proof, and some care should be executed when large basis 
sets are being used. 

As mentioned above, it is good policy to treat the primitive 
functions in groups, whether they are to be contracted or not. 

The number of final , contracted functions in such a group should 
not be too large. For instance, as an example, a contraction 
group involving e.g. three contracted functions of F-type 
(irrespective of the number of primitives) would require 3 4 *10 4 = 
810 000 words of storage in order to accomodate the contracted 
f f f f-integrals . treat the primitive functions in groups, whether 
they are to be contracted or not. The number of final , contracted 
functions in such a group should not be too large. For instance, 
as an example, a contraction group involving e.g. three 
contracted functions of F-type (irrespective of the number of 
primitives) would require 3 4 »10 4 = 810 000 words of storage in 
order to accomodate the contracted fff f-integrals . 


In most cases, the alarm for insufficient memory will not be 
alerted until an integral block which violates memory limitations 
is treated. For obvious reasons, it is wise to place the heavy 
atoms first in the input sequence, such that a possible crash is 
encountered early in the calculation. 


9. 


Plans for future development 


The present code is supposed to be a first prototype useful for 
production calculations. A number of changes and additions are 
planned, however, which are going to increase its usefulness in 
different respects. Below, we list a number of the features which 
are under development: 

a) Special code for low quantum numbers. See Section 5 or a 
further discussion on this topic. 

b) New strategy for high quantum numbers with few primitives. 
The present code performs best for angular quantum numbers in the 
medium range (s, p, and d functions), and when the number of 
primitives is sufficient to allow for efficient vectorization . 

For solitary primitive functions with high L there are better 
strategies, which vectorize over sums of L-values rather than 
over primitive integrals. We are presently working on a code 
using the method of McMurchie and Davidson /5/, which will be 
addeed to the program in the future. 

c) Vectorization over symmetry- related integrals. 

In many cases it will be possible to increase the effective vec- 
tor length - and thus the performance of the program - by treat- 
ing some or all symmetry-related integrals in the same batch. 

This will require relatively few changes of the code, and steps 
are already being taken to prepare for such a modification. 


d) Direct SCF implementation. 

The Direct SCF technique /6/ offers interesting possibilities for 
very large SCF calculations, particularly in a combination with 
the CPU-power available on modern hardware. However, the present 
version, as implemented in the program DISCO, uses the old MOLE- 
CULE integral routines. Replacing these by the newly written code 
will drastically increase the size of molecules that can be 
studied with this technique. 
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Appendix A. Definitions and abbreviations. 


1) A GTO (Gaussian Type Orbital) is a primitive Cartesian Gaus- 
sian function: 

X A (r) = (x-x A ) 1 (y-y A )^(z-z A ) k exp(-A(r-r A ) 2 ) (1) 

2) A CGTO (Contracted GTO) is a linear combination of GTOs on 
the same center r A and with the same 'quantum' numbers i, j, 
and k (i.e. differing in the orbital exponent A only). 

3) A GTS is a shell of GTOs, i.e. a set of GTOs differing only 

in the quantum numbers i, j and k. A complete set is made up 

by all combinations of i, j and k subject to the restriction 

i+j+k = L. This leads to a set of < L J 2 ) functions, or 1, 3, 
6, 10, and 15 functions for L = 0, 1, 2, 3, or 4. 

4) A CGTS is a similar shell of CGTOs. 

5) A GGTO is a group of GTOs used in a contraction as described 

in 2 ) above . 

6) A GGTS , similarly, is a group of GTSs (i.e. a shell of 
GGTOs ) . 



Appendix B. Input description for MOLECULE (1984 version). 
All input is given in free format. 


Record 1 : 


Record 2 : 


Arbitrary title. 


NTYP , NSOP , ( SYM ( I ) , 1=1 ,NSOP ) , XLIM, IUNIT, IPR 


NTYP = number of different types of atoms (i.e. 
number of different basis sets). 

NSOP = number of symmetry elements generating the 
molecular point group to be used. 

SYM (type CHARACTER) defines the symmetry elements. 

XLIM = threshold for the integral calculation. 

IUNIT = 0 assumes geometry input in atomic units. 

” =1 " *’ " Angstrom. 

IPR = print parameter. IPR=0 gives a minimum of 
printout, IPR=3 gives maximum. 



Group 3 : 


To be repeated NTYP tines. 


Record 3 . 


Note: 


: BAS , NON , Q , SPDF , ( JCO ( J ) , J=1 , 5 ) 

BAS (type CHARACTER) = identifier of a standard 
basis set to be read fron a library (this will need 
some local modifications to work). 

If BAS = 'blank' the basis set is read from the 
input file. 

NON = number of symmetry-independent atoms using 
this basis set. 

Q = nuclear charge. 

SPDF (type CHARACTER) defines the highest L-value in 
this basis set. 

JCO = number of groups of primitive GTOs for each 
L-value . 

Primitive GTOs in the same contraction are always in 
the same group. For efficiency reasons, one may also 
group together primitives that are not to be con- 
tracted. The performance of the code is related to 
the degree of grouping. Single GTOs should be avoid- 
ed whenever possible. 



Record 3.2: To be repeated NON times. 


NAME , X , Y , Z 

NAME (type CHARACTER) = label for each symmetry- 
independent atom. 

X, Y, Z = Cartesian coordinates for each symmetry- 
independent atom. 

Group 3.3: To be repeated for each L- value. 

(Group 3.3 is to be specified only if BAS = 
'blank' in Record 3.1). 

Record 3.3.1: ( NUC(K) , K=1,JC0(J)) 

NUC (K) = number of primitive GTOs in group K. 

Record 3.3.2: ( NRC(K) , K=1,JC0(J)) 

NRC (K) = number of contracted functions (CGTOs) 


in group K. 



Group 3.3.3: 


To be repeated JCO(J) times. 


Record 3. 3. 3.1: (ALPHA (N ) , N=1,NUC(K)) 

orbital exponents for the contraction group K. 


Record 3. 3. 3. 2: To be repeated NRC(K) times. 


(CONT(N), N=1 , NUC ( K ) ) 


contraction coefficients for the group K. 


